# Directory and data import
setwd("/home/joshk/git_repositories/BrookTrout_TG/Data")
Trout <- read.csv("BrookTrout.csv")
Trout <- subset(Trout, !is.na(Mother.ID))
Trout <- subset(Trout, !is.na(Father.ID))
Trout$Offspring.Acc <- as.factor(Trout$Offspring.Acc)
Trout$Mother.Acc <- as.factor(Trout$Mother.Acc)
Trout$Father.Acc <- as.factor(Trout$Father.Acc)
Trout$Mother.ID <- as.factor(Trout$Mother.ID)
Trout$Father.ID <- as.factor(Trout$Father.ID)
options(na.action = "na.fail")
## Loading in necessary packages
library("easypackages")
libraries(
"car", "dotwhisker", "dplyr", "ggplot2", "grid", "gridExtra",
"knitr", "lme4", "magrittr", "MuMIn", "nlme", "purrr", "rmarkdown"
)Removing mass-specific MO2 measurements from data-frame
Data = Trout[,c(1:10, 14:17)]Assigning ggplot2 themes for figures
my.theme <- theme(
panel.grid.minor = element_blank(),
axis.title = element_text(size = 14, family = "Noto Sans"),
axis.text = element_text(size = 14, colour = "black", family = "Noto Sans"),
axis.title.y = element_text(angle = 0, vjust = 0.5), panel.grid.major =
element_line(colour = "grey75"), legend.title = element_text(
size = 14,
colour = "black", family = "Noto Sans"
), legend.text = element_text(
size = 12,
colour = "black", family = "Noto Sans"
)
)
my.theme.dw <- theme(
panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
axis.title = element_text(size = 12, family = "Noto Sans"),
axis.text = element_text(size = 12, colour = "black", family = "Noto Sans"),
legend.title = element_text(
size = 12,
colour = "black", family = "Noto Sans"
), legend.text = element_text(
size = 12,
colour = "black", family = "Noto Sans"
)
)
my.theme.diag <- theme(
panel.grid.minor = element_blank(), panel.grid.major = element_blank(),
axis.title = element_text(size = 8, family = "Noto Sans"),
axis.text = element_text(size = 8, colour = "black", family = "Noto Sans"),
legend.title = element_text(
size = 8,
colour = "black", family = "Noto Sans"
), legend.text = element_text(
size = 8,
colour = "black", family = "Noto Sans"
)
)First, assessing distribution of response variables
{
RMR.Hist <- ggplot(data = Data, aes(x = RMR)) +
geom_histogram(colour = "black", fill = "royalblue") +
my.theme +
theme_bw() +
geom_vline(
xintercept = with(na.omit(Data), mean(RMR)), size = 1.25,
colour = "darkred", linetype = "dotted"
) +
xlab("Resting Metabolic \n Rate (mg O2/h)")
MMR.Hist <- ggplot(data = Data, aes(x = MMR)) +
geom_histogram(colour = "black", fill = "royalblue") +
my.theme +
theme_bw() +
geom_vline(
xintercept = with(na.omit(Data), mean(MMR)), size = 1.25,
colour = "darkred", linetype = "dotted"
) +
xlab("Peak Metabolic \n Rate (mg O2/h)")
Mass.hist <- ggplot(data = Trout, aes(x = Mass)) +
geom_histogram(colour = "black", fill = "royalblue") +
my.theme +
theme_bw() +
geom_vline(
xintercept = with(Trout, mean(Mass)), size = 1.25,
colour = "darkred", linetype = "dotted"
) +
xlab("Mass")
CTM.hist <- ggplot(data = Trout, aes(x = CTM)) +
geom_histogram(colour = "black", fill = "royalblue") +
my.theme +
theme_bw() +
geom_vline(
xintercept = with(subset(Data, !is.na(CTM)), mean(CTM)), size = 1.25,
colour = "darkred", linetype = "dotted"
) +
xlab("Critical Thermal Maximum (Degrees C)")
}
grid.arrange(RMR.Hist, MMR.Hist, Mass.hist, CTM.hist, nrow = 2)
## Mass appears to fit a Gaussian distibution. Critical thermal maximum not awfully skewed, and all non-mass specific parameters meet gaussian distributions.Checking for patterns or biases in mass by family and acclimation temperature
# Boxplots
{
Mass.plot1 <- ggplot(data = Data, aes(x = as.factor(Parent.group), y = Mass)) +
my.theme +
theme_bw() +
geom_boxplot(fill = "cadetblue") +
xlab("Parent Group")
Mass.plot2 <- ggplot(data = Data, aes(x = as.factor(Family), y = Mass)) +
my.theme +
theme_bw() +
geom_boxplot(fill = "mediumorchid") +
xlab("Family")
Maternal.Bin <- ggplot(data = Data, aes(x = as.factor(Mother.Acc), y = Mass)) +
my.theme +
theme_bw() +
geom_boxplot(fill = "gold2") +
xlab("Maternal Acclimation")
Paternal.Bin <- ggplot(data = Data, aes(x = as.factor(Father.Acc), y = Mass)) +
my.theme +
theme_bw() +
geom_boxplot(fill = "darkseagreen2") +
xlab("Paternal Acclimation")
Offspring.Bin <- ggplot(data = Data, aes(x = as.factor(Offspring.Acc), y = Mass)) +
my.theme +
theme_bw() +
geom_boxplot(fill = "lightcoral") +
xlab("Offspring Acclimation")
Details <- ggplot(data = data.frame("x" = c(1:100), "y" = c(1:100)), aes(x = x, y = y)) +
theme_void() +
theme(legend.position = "none") +
annotate("text",
x = 25, y = 70, label =
"Central lines represent the \n median of each category."
)
}
grid.arrange(Mass.plot1, Mass.plot2, Maternal.Bin, Paternal.Bin, Offspring.Bin,
Details,
nrow = 3, top = "Mass by Group"
)
## Large variation between families, and significant variation between parent groups. Surprising
## that offspring from cold acclimated parents are larger, but not surprising that warm acclimated
## offspring are larger.
{
RMR.by.Treat <- ggplot(data = Data, aes(x = as.factor(Parent.group), y = RMR)) +
my.theme +
theme_bw() +
geom_boxplot(fill = "cadetblue") +
xlab("Parent Group") +
ylab("Resting Metabolic Rate \n (mg O2/h)")
MMR.by.Treat <- ggplot(data = Data, aes(x = as.factor(Parent.group), y = AAS)) +
my.theme +
theme_bw() +
geom_boxplot(fill = "palegreen3") +
xlab("Parent Group") +
ylab("Peak Metabolic Rate \n (mg O2/h)")
CTM.by.Treat <- ggplot(data = Data, aes(x = as.factor(Parent.group), y = CTM)) +
my.theme +
theme_bw() +
geom_boxplot(fill = "thistle") +
xlab("Parent Group") +
ylab("Critical Thermal Maximum \n (Degrees C)")
}
grid.arrange(RMR.by.Treat, MMR.by.Treat, CTM.by.Treat,
nrow = 2, top = "Metabolic Rate and Critical Thermal Max by Treatment"
)
# No distinct trends outside of peak metabolic rate by parent groups 1 and 4.
{
RMR.by.Fam <- ggplot(data = Data, aes(x = as.factor(Family), y = RMR)) +
my.theme +
theme_bw() +
geom_boxplot(fill = "cadetblue") +
xlab("Parent Group") +
ylab("Resting Metabolic Rate \n (mg O2/h)")
MMR.by.Fam <- ggplot(data = Data, aes(x = as.factor(Family), y = MMR)) +
my.theme +
theme_bw() +
geom_boxplot(fill = "palegreen3") +
xlab("Family") +
ylab("Peak Metabolic Rate \n (mg O2/h)")
CTM.by.Fam <- ggplot(data = Data, aes(x = as.factor(Family), y = CTM)) +
my.theme +
theme_bw() +
geom_boxplot(fill = "thistle") +
xlab("Family") +
ylab("Critical Thermal Maximum \n (Degrees C)")
}
grid.arrange(RMR.by.Fam, MMR.by.Fam, CTM.by.Fam,
nrow = 2, top = "Metabolic Rate and Critical Thermal Max by Family"
)
## No consistent patterns accross families.Exploring effects of mass on metabolic parameters with preliminary plots
{
RMR.by.Mass <- ggplot(data = Data, aes(x = Mass, y = RMR)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "cadetblue", alpha = 0.3) +
xlab("Mass (g)") +
geom_smooth(method = lm, colour = "black") +
ylab("Resting Metabolic Rate \n (mg O2/h)")
MMR.by.Mass <- ggplot(data = Data, aes(x = Mass, y = MMR)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "thistle", alpha = 0.3) +
xlab("Mass (g)") +
geom_smooth(method = lm, colour = "black") +
ylab("Peak Metabolic Rate \n (mg O2/h)")
CTM.by.Mass <- ggplot(data = Data, aes(x = Mass, y = CTM)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "gold2", alpha = 0.3) +
xlab("Mass (g)") +
geom_smooth(method = lm, colour = "black") +
ylab("Critical Thermal Maximum \n (Degrees C)")
}
grid.arrange(RMR.by.Mass, MMR.by.Mass, CTM.by.Mass,
nrow = 2, top = "Metabolic Rate and Critical Thermal Maximum by Mass"
)
# No clear outliers, but critical thermal maximum and resting metabolic rate heterogenous
# by mass. Weighting may be necessary if mass-specific responses not used.Checking for outliers with Cleveland dotplots
{
RMR.NMS.Data <- subset(Data, !is.na(RMR))
RMR.dotplot <- ggplot(data = RMR.NMS.Data, aes(
x = RMR, y =
(1:nrow(RMR.NMS.Data))
)) +
my.theme +
theme_bw() +
geom_point(
size =
3, colour = "gold2", alpha = 0.5
) +
xlab("Resting Metabolic Rate (mg O2/h)") +
ylab("Row Number") +
annotate("rect",
xmin = (with(RMR.NMS.Data, mean(RMR) - 3 * sd(RMR))),
xmax = (with(RMR.NMS.Data, mean(RMR) + 3 * sd(RMR))), ymin = 0, ymax = nrow(RMR.NMS.Data),
alpha = .2
)
MMR.NMS.Data <- subset(Data, !is.na(MMR))
MMR.dotplot <- ggplot(data = MMR.NMS.Data, aes(
x = MMR, y =
(1:nrow(MMR.NMS.Data))
)) +
my.theme +
theme_bw() +
geom_point(
size =
3, colour = "royalblue", alpha = 0.5
) +
xlab("Peak Metabolic Rate (mg O2/h)") +
ylab("Row Number") +
annotate("rect",
xmin = (with(MMR.NMS.Data, mean(MMR) - 3 * sd(MMR))),
xmax = (with(MMR.NMS.Data, mean(MMR) + 3 * sd(MMR))), ymin = 0, ymax = nrow(MMR.NMS.Data),
alpha = .2
)
CTM.NMS.Data <- subset(Data, !is.na(CTM))
CTM.dotplot <- ggplot(data = CTM.NMS.Data, aes(
x = CTM, y =
(1:nrow(CTM.NMS.Data))
)) +
my.theme +
theme_bw() +
geom_point(
size =
3, colour = "black", alpha = 0.5
) +
xlab("Critical Thermal Maximum (Degrees C)") +
ylab("Row Number") +
annotate("rect",
xmin = (with(CTM.NMS.Data, mean(CTM) - 3 * sd(CTM))),
xmax = (with(CTM.NMS.Data, mean(CTM) + 3 * sd(CTM))), ymin = 0, ymax = nrow(CTM.NMS.Data),
alpha = .2
)
}
grid.arrange(RMR.dotplot, MMR.dotplot, CTM.dotplot,
nrow = 2, top = "Non-Mass Specific Response Dotplots: \n Grey boxed represent mean +/-
three standard deviations"
)
## Outliers not extremely distinct, so retained but monitored below.Non-Mass Specific Linear Models
Note that sires are not unique to dams; therefore random effects should be crossed
{
RMR.model <- lmer(RMR ~ Mother.Acc * Father.Acc * Offspring.Acc + Mass + (1 | Mother.ID) + (1 | Father.ID),
REML = FALSE, data = RMR.NMS.Data
)
MMR.model <- lmer(MMR ~ Mother.Acc * Father.Acc * Offspring.Acc + Mass + (1 | Mother.ID) + (1 | Father.ID),
REML = FALSE, data = MMR.NMS.Data
)
CTM.model <- lmer(CTM ~ Mother.Acc * Father.Acc * Offspring.Acc + Mass + (1 | Mother.ID) + (1 | Father.ID),
REML = FALSE, weights = Mass, data = CTM.NMS.Data
)
}
# Checking sample sizes
paste("Resting MO2", nrow(RMR.NMS.Data), sep = ": ")## [1] "Resting MO2: 186"paste("Peak MO2", nrow(MMR.NMS.Data), sep = ": ")## [1] "Peak MO2: 224"paste("Critical Thermal Max", nrow(CTM.NMS.Data), sep = ": ")## [1] "Critical Thermal Max: 216"Quick assessment of model dispersion
{
RMR.Dev <- abs(summary(RMR.model)$AICtab[[4]] / summary(RMR.model)$AICtab[[5]])
MMR.Dev <- abs(summary(MMR.model)$AICtab[[4]] / summary(MMR.model)$AICtab[[5]])
CTM.Dev <- abs(summary(CTM.model)$AICtab[[4]] / summary(CTM.model)$AICtab[[5]])
}
print(paste("RMR:", RMR.Dev, "MMR:", MMR.Dev, "CTM:", CTM.Dev, sep = " "))## [1] "RMR: 0.363211108441395 MMR: 0.0141738136619573 CTM: 0.0593597551129786"# All models are highly over-dispersedRepeating models with log-transformed (here, using the natural log) response variable and mass
{
RMR.model <- lmer(log(RMR) ~ Mother.Acc * Father.Acc * Offspring.Acc + log(Mass) + (1 | Mother.ID) + (1 | Father.ID),
REML = FALSE, data = RMR.NMS.Data
)
MMR.model <- lmer(log(MMR) ~ Mother.Acc * Father.Acc * Offspring.Acc + log(Mass) + (1 | Mother.ID) + (1 | Father.ID),
REML = FALSE, data = MMR.NMS.Data
)
CTM.model <- lmer(CTM ~ Mother.Acc * Father.Acc * Offspring.Acc + log(Mass) + (1 | Mother.ID) + (1 | Father.ID),
REML = FALSE, weights = Mass, data = CTM.NMS.Data
)
}
# Double-checking dispersion
{
RMR.Dev <- abs(summary(RMR.model)$AICtab[[4]] / summary(RMR.model)$AICtab[[5]])
MMR.Dev <- abs(summary(MMR.model)$AICtab[[4]] / summary(MMR.model)$AICtab[[5]])
CTM.Dev <- abs(summary(CTM.model)$AICtab[[4]] / summary(CTM.model)$AICtab[[5]])
}
print(paste("RMR:", RMR.Dev, "MMR:", MMR.Dev, "CTM:", CTM.Dev, sep = " "))## [1] "RMR: 1.29944000421038 MMR: 0.045629817986151 CTM: 0.0595168880245248"## Dispersion slightly better, particularly for resting MO2.
# Checking sample size for each grouping, and calculating statistical power of models
RMR.NMS.Data %>%
group_by(Offspring.Acc, Father.Acc, Mother.Acc) %>%
summarise(N = n()) %>%
as.data.frame() %>%
summarise(Mean = mean(N), SD = sd(N))## Mean SD
## 1 23.25 5.700877require("pwr")
num_df <- nrow(summary(RMR.model)$coefficients) - 1
den_df <- nrow(RMR.NMS.Data) - (num_df + 1)
pwr.f2.test(u = num_df, v = den_df, f2 = 0.05, sig.level = 0.05)##
## Multiple regression power calculation
##
## u = 8
## v = 177
## f2 = 0.05
## sig.level = 0.05
## power = 0.5284571# Fairly high statistical power here. Iterating models and sorting final model outputs by AICc
AIC.Results.NMS <- vector("list", 3)
models <- c("RMR.model", "MMR.model", "CTM.model")
# Looping through models, running with all combinations of predictors,
# then sorting by AICc
for (i in 1:length(models)) {
AIC.Results.NMS[[i]] <- subset(dredge(get(models[i]),
rank = "AICc",
extra = list(R2 = function(x) r.squaredGLMM(x))
), delta < 5)
}
Select.Models.NMS <- c()
Delta.NMS <- c()
Response <- c()
Marg.R2 <- c()
Cond.R2 <- c()
# Looping through models and collecting metrics of fit (i.e. marginal and conditional r-squared values)
for (j in 1:length(AIC.Results.NMS)) {
for (i in 1:length(attr(AIC.Results.NMS[[j]], "model.calls"))) {
Select.Models.NMS[i] <- gsub(".*~ ", "", paste(attr(AIC.Results.NMS[[j]], "model.calls")[[i]])[2])
Delta.NMS[i] <- AIC.Results.NMS[[j]]$delta[i]
Response[i] <- gsub(" ~.*", "", paste(attr(AIC.Results.NMS[[j]], "model.calls")[[i]])[2])
Marg.R2[i] <- AIC.Results.NMS[[j]]$R21[i]
Cond.R2[i] <- AIC.Results.NMS[[j]]$R22[i]
if (i == length(attr(AIC.Results.NMS[[j]], "model.calls"))) {
Table <- data.frame(
"Reponse.Type" = "Mass Independant", "Response" = Response,
"Delta.AIC" = Delta.NMS, "Marginal.R2" = Marg.R2, "Conditional.R2" = Cond.R2,
"Equation" = Select.Models.NMS
)
assign(paste("AIC.Selection", "NMS", j, sep = "_"), Table)
}
}
Select.Models.NMS <- c()
Delta.NMS <- c()
Response <- c()
Marg.R2 <- c()
Cond.R2 <- c()
}
# Binding outcomes together
AIC.Selection.NMS <- rbind(AIC.Selection_NMS_1, AIC.Selection_NMS_2, AIC.Selection_NMS_3)
# Viewing and saving output
View(AIC.Selection.NMS)
write.csv(AIC.Selection.NMS, "Model_AIC_Final.csv", row.names = FALSE)Assessing model outputs
Testing model assumptions for criticial thermal maximum: Model 1.
## First assessing response by predictors;
{
PatT.CTM <- ggplot(data = CTM.NMS.Data, aes(
x = as.factor(Father.Acc),
y = CTM
)) +
theme_bw() +
geom_boxplot(fill = "mediumorchid", aes(middle = mean(CTM))) +
xlab("Sire Acclimation Temperature") +
ylab("Critical Thermal Maximum \n (Degrees C)")
DamT.CTM <- ggplot(data = CTM.NMS.Data, aes(
x = as.factor(Mother.Acc),
y = CTM
)) +
theme_bw() +
geom_boxplot(fill = "mediumseagreen", aes(middle = mean(CTM))) +
xlab("Dam Acclimation Temperature") +
ylab("Critical Thermal Maximum \n (Degrees C)")
OffT.CTM <- ggplot(data = CTM.NMS.Data, aes(
x = as.factor(Offspring.Acc),
y = CTM
)) +
theme_bw() +
geom_boxplot(fill = "royalblue", aes(middle = mean(CTM))) +
xlab("Offspring Acclimation Temperature") +
ylab("Critical Thermal Maximum \n (Degrees C)")
Pat.by.CTM1 <- ggplot(data = CTM.NMS.Data, aes(
x = as.factor(Father.Acc),
y = CTM, fill = Offspring.Acc
)) +
theme_bw() +
theme(legend.position = "top") +
geom_boxplot() +
scale_fill_manual(values = c("mediumorchid", "royalblue")) +
xlab("Sire Acclimation Temperature") +
ylab("Critical Thermal Maximum \n (Degrees C)")
Pat.by.CTM2 <- ggplot(data = CTM.NMS.Data, aes(
x = as.factor(Father.Acc),
y = CTM, fill = Mother.Acc
)) +
theme_bw() +
theme(legend.position = "top") +
geom_boxplot() +
scale_fill_manual(values = c(
"mediumseagreen",
"royalblue"
)) +
xlab("Sire Acclimation Temperature") +
ylab("Critical Thermal Maximum \n (Degrees C)")
Off.by.CTM <- ggplot(data = CTM.NMS.Data, aes(
x = as.factor(Offspring.Acc),
y = CTM, fill = Mother.Acc
)) +
theme_bw() +
theme(legend.position = "top") +
geom_boxplot() +
scale_fill_manual(values = c(
"mediumorchid",
"mediumseagreen"
)) +
xlab("Offspring Acclimation Temperature") +
ylab("Critical Thermal Maximum \n (Degrees C)")
}
grid.arrange(PatT.CTM, DamT.CTM, OffT.CTM, Pat.by.CTM1, Pat.by.CTM2, Off.by.CTM, nrow = 2, top = "Fixed Effects by Response: Boxplot Bars Represent Means")
## Re-iterating top model for further plotting
CTM.m1 <- lmer(CTM ~ Father.Acc + Mother.Acc + Offspring.Acc + Father.Acc:Mother.Acc +
Father.Acc:Offspring.Acc + Mother.Acc:Offspring.Acc +
Father.Acc:Mother.Acc:Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID),
REML = FALSE, weights = Mass, data = CTM.NMS.Data
)
# Note singular fit error. Checking estimates of random intercepts
summary(CTM.m1)$varcor## Groups Name Std.Dev.
## Mother.ID (Intercept) 2.0348e-02
## Father.ID (Intercept) 9.2761e-09
## Residual 4.0160e-01## Paternal ID explaining little to no variance in data. Could probably be dropped from model. First, diagnosing residuals.
{
Res.Alone.m1 <- ggplot(
data = data.frame("Res" = residuals(CTM.m1, type = "deviance")),
aes(x = 1:nrow(CTM.NMS.Data), y = Res)
) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) +
xlab("Row Number") +
ylab("Deviance Residuals")
Res.by.Fit.m1 <- ggplot(
data = data.frame("Res" = residuals(CTM.m1), "Fit" = predict(CTM.m1)),
aes(x = Fit, y = Res)
) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) +
xlab("Fitted Values") +
ylab("Deviance Residuals")
Res.by.Resp.m1 <- ggplot(
data = data.frame("Res" = residuals(CTM.m1), "Resp" = CTM.NMS.Data$CTM),
aes(x = Resp, y = Res)
) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) +
xlab("log Mass Specific Metabolic Rate") +
ylab("Deviance Residuals")
QQNorm.1.m1 <- ggplot(
data.frame("Res" = residuals(CTM.m1), "Fit" = predict(CTM.m1)),
aes(sample = Res)
) +
stat_qq(colour = "gold") +
stat_qq_line() +
my.theme +
theme_bw()
}
grid.arrange(Res.Alone.m1, Res.by.Fit.m1, Res.by.Resp.m1, QQNorm.1.m1, nrow = 2, top = "Mass Specific Resting Metabolic Rate Residuals")
qqPlot(residuals(CTM.m1, type = "deviance"), ylab = "Deviance Residuals")
## 102 14
## 97 14# Some low end deviation, but minimal.
CTM.m1.res <- data.frame("Model" = "Model1", "Res" = c(residuals(CTM.m1, type = "deviance")))
ggplot(CTM.m1.res, aes(Res, fill = Model, colour = Model)) +
geom_histogram(alpha = 0.5) +
theme_bw() +
my.theme +
xlab("Deviance Residuals") +
ylab("Count") +
scale_fill_manual(
values =
c("mediumseagreen")
) +
scale_colour_manual(values = c("black")) +
theme(legend.position = "none") +
geom_density(alpha = 0.5, aes(y = 0.1 * ..count..))
# Residuals appear reasonable.
## Assessing parameter-specific patterns.
{
PatT.Res.m1 <- ggplot(data = CTM.NMS.Data, aes(
x = as.factor(Father.Acc),
y = residuals(CTM.m1, "deviance")
)) +
theme_bw() +
geom_boxplot(fill = "mediumorchid", aes(middle = mean(residuals(CTM.m1, "deviance")))) +
xlab("Sire Acclimation Temperature") +
ylab("Deviance Residuals")
DamT.Res.m1 <- ggplot(data = CTM.NMS.Data, aes(
x = as.factor(Mother.Acc),
y = residuals(CTM.m1, "deviance")
)) +
theme_bw() +
geom_boxplot(fill = "mediumseagreen", aes(middle = mean(residuals(CTM.m1, "deviance")))) +
xlab("Dam Acclimation Temperature") +
ylab("Deviance Residuals")
OffT.Res.m1 <- ggplot(data = CTM.NMS.Data, aes(
x = as.factor(Offspring.Acc),
y = residuals(CTM.m1, "deviance")
)) +
theme_bw() +
geom_boxplot(fill = "royalblue", aes(middle = mean(residuals(CTM.m1, "deviance")))) +
xlab("Offspring Acclimation Temperature") +
ylab("Deviance Residuals")
Pat.by.Off.T.Res.m1 <- ggplot(data = CTM.NMS.Data, aes(
x = as.factor(Father.Acc),
y = residuals(CTM.m1, type = "deviance"), fill = Offspring.Acc
)) +
theme_bw() +
theme(legend.position = "top") +
geom_boxplot() +
scale_fill_manual(values = c("mediumorchid", "royalblue")) +
xlab("Sire Acclimation Temperature") +
ylab("Deviance Residuals")
Pat.by.Dam.T.Res.m1 <- ggplot(data = CTM.NMS.Data, aes(
x = as.factor(Father.Acc),
y = residuals(CTM.m1, type = "deviance"), fill = Mother.Acc
)) +
theme_bw() +
theme(legend.position = "top") +
geom_boxplot() +
scale_fill_manual(values = c(
"mediumseagreen",
"royalblue"
)) +
xlab("Sire Acclimation Temperature") +
ylab("Deviance Residuals")
Off.by.Dam.T.Res.m1 <- ggplot(data = CTM.NMS.Data, aes(
x = as.factor(Offspring.Acc),
y = residuals(CTM.m1, type = "deviance"), fill = Mother.Acc
)) +
theme_bw() +
theme(legend.position = "top") +
geom_boxplot() +
scale_fill_manual(values = c(
"mediumorchid",
"mediumseagreen"
)) +
xlab("Offspring Acclimation Temperature") +
ylab("Deviance Residuals")
}
grid.arrange(PatT.Res.m1, DamT.Res.m1, OffT.Res.m1, Pat.by.Off.T.Res.m1, Pat.by.Dam.T.Res.m1, Off.by.Dam.T.Res.m1, nrow = 2, top = "Fixed Effects by Residuals: Boxplot Bars Represent Means")
# Plots suggest a possible sire by offspring heteroskedasticity. Retrying model with constant variance per group.
CTM.m1.1 <- lme(CTM ~ Father.Acc + Mother.Acc + Offspring.Acc + Father.Acc:Mother.Acc +
Father.Acc:Offspring.Acc + Mother.Acc:Offspring.Acc +
Father.Acc:Mother.Acc:Offspring.Acc,
random = list(~ 1 | Mother.ID, ~ 1 | Father.ID),
weights = ~Mass, data = CTM.NMS.Data
)
CTM.m1.2 <- lme(CTM ~ Father.Acc + Mother.Acc + Offspring.Acc + Father.Acc:Mother.Acc +
Father.Acc:Offspring.Acc + Mother.Acc:Offspring.Acc +
Father.Acc:Mother.Acc:Offspring.Acc,
random = list(~ 1 | Mother.ID, ~ 1 | Father.ID),
weights = varIdent(form = ~ Mass | Father.Acc * Offspring.Acc), data = CTM.NMS.Data
)
# Note that singular fit error is no longer present.
summary(CTM.m1.1)$AIC / summary(CTM.m1.2)$AIC## [1] 6.157988# Wow! Significant improvement in AIC. Again, assessing residuals.
{
Res.Alone.m1 <- ggplot(
data = data.frame("Res" = residuals(CTM.m1.2, type = "normalized")),
aes(x = 1:nrow(CTM.NMS.Data), y = Res)
) +
theme_bw() +
geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) +
xlab("Row Number") +
ylab("Deviance Residuals")
Res.by.Fit.m1 <- ggplot(data = data.frame(
"Res" = residuals(CTM.m1.2, type = "normalized"),
"Fit" = predict(CTM.m1.2)
), aes(x = Fit, y = Res)) +
theme_bw() +
geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) +
xlab("Fitted Values") +
ylab("Deviance Residuals")
Res.by.Resp.m1 <- ggplot(data = data.frame(
"Res" = residuals(CTM.m1.2, type = "normalized"),
"Resp" = CTM.NMS.Data$CTM
), aes(x = Resp, y = Res)) +
theme_bw() +
geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) +
xlab("log Mass Specific Metabolic Rate") +
ylab("Deviance Residuals")
QQNorm.1.m1 <- ggplot(data.frame(
"Res" = residuals(CTM.m1.2, type = "normalized"),
"Fit" = predict(CTM.m1.2)
), aes(sample = Res)) +
stat_qq(colour = "gold") +
stat_qq_line() +
theme_bw()
}
grid.arrange(Res.Alone.m1, Res.by.Fit.m1, Res.by.Resp.m1, QQNorm.1.m1, nrow = 2, top = "Mass Specific Resting Metabolic Rate Residuals")
# Far cleaner. Double-checking confidence intervals.
qqPlot(residuals(CTM.m1.2, type = "normalized"), ylab = "Normalized Residuals")
## 1/3 1/3
## 154 161# Well enough behaved! Re-testing AIC with new variance structure.
CTM.model <- lme(CTM ~ Mother.Acc * Father.Acc * Offspring.Acc + Mass,
random = list(~ 1 | Mother.ID, ~ 1 | Father.ID),
weights = varIdent(~ Mass | Father.Acc * Offspring.Acc), data = CTM.NMS.Data
)
subset(dredge(CTM.model, rank = "AICc"), delta <= 2.0)## Global model call: lme.formula(fixed = CTM ~ Mother.Acc * Father.Acc * Offspring.Acc +
## Mass, data = CTM.NMS.Data, random = list(~1 | Mother.ID,
## ~1 | Father.ID), weights = varIdent(~Mass | Father.Acc *
## Offspring.Acc))
## ---
## Model selection table
## (Int) Off.Acc df logLik AICc delta weight
## 9 28.57 + 5 -19.117 48.5 0 1
## Models ranked by AICc(x)
## Random terms (all models):
## '1 | Mother.ID', '1 | Father.ID %in% Mother.ID'Re-calculating AICs for model iterations: criticial thermal maximum.
AIC.Results.CTM1 <- subset(dredge(CTM.model,
rank = "AICc",
extra = list(R2 = function(x) r.squaredGLMM(x))
), delta < 2)
{
Select.Models.CTM <- c()
Delta.NMS <- c()
Response <- c()
Marg.R2 <- c()
Cond.R2 <- c()
Select.Models.CTM1 <- gsub("[[:space:]].*", "", gsub(".*~ ", "", paste(attr(AIC.Results.CTM1, "model.calls")[[1]])[2]))
Delta.CTM1 <- AIC.Results.CTM1$delta[1]
Response <- gsub(" ~.*", "", paste(attr(AIC.Results.CTM1, "model.calls")[[1]])[2])
Marg.R2 <- AIC.Results.CTM1$R21[1]
Cond.R2 <- AIC.Results.CTM1$R22[1]
Table <- data.frame("Reponse.Type" = "Mass Independant", "Response" = Response, "Delta.AIC" = Delta.CTM1, "Marginal.R2" = Marg.R2, "Conditional.R2" = Cond.R2, "Equation" = Select.Models.CTM1)
assign(paste("AIC.Selection", "NMS", 3, sep = "_"), Table)
}
# Stacking all top models
AIC.Selection.NMS <- rbind(AIC.Selection_NMS_1, AIC.Selection_NMS_2, AIC.Selection_NMS_3)
View(AIC.Selection.NMS)
write.csv(AIC.Selection.NMS, "Model_AIC_Final.csv", row.names = FALSE)Re-assessing criticial thermal maximum: model 1.
CTM.model.m1 <- lme(CTM ~ Offspring.Acc, random = list(~ 1 | Mother.ID, ~ 1 | Father.ID), weights = varIdent(~ Mass | Father.Acc * Offspring.Acc), data = CTM.NMS.Data)
# Residual Plots
{
Res.Alone.m1 <- ggplot(
data = data.frame("Res" = residuals(CTM.model.m1, type = "normalized")),
aes(x = 1:nrow(CTM.NMS.Data), y = Res)
) +
theme_bw() +
geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) +
xlab("Row Number") +
ylab("Deviance Residuals")
Res.by.Fit.m1 <- ggplot(data = data.frame(
"Res" = residuals(CTM.model.m1, type = "normalized"),
"Fit" = predict(CTM.model.m1)
), aes(x = Fit, y = Res)) +
theme_bw() +
geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) +
xlab("Fitted Values") +
ylab("Deviance Residuals")
Res.by.Resp.m1 <- ggplot(data = data.frame(
"Res" = residuals(CTM.model.m1, type = "normalized"),
"Resp" = CTM.NMS.Data$CTM
), aes(x = Resp, y = Res)) +
theme_bw() +
geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) +
xlab("log Mass Specific Metabolic Rate") +
ylab("Deviance Residuals")
QQNorm.1.m1 <- ggplot(data.frame(
"Res" = residuals(CTM.model.m1, type = "normalized"),
"Fit" = predict(CTM.model.m1)
), aes(sample = Res)) +
stat_qq(colour = "gold") +
stat_qq_line() +
theme_bw()
}
grid.arrange(Res.Alone.m1, Res.by.Fit.m1, Res.by.Resp.m1, QQNorm.1.m1,
nrow = 2, top =
"Mass Specific Resting Metabolic Rate Residuals"
)
qqPlot(residuals(CTM.model.m1, type = "normalized"), ylab = "Deviance Residuals")
## 1/1 3/3
## 14 89# Checking marginal and conditional R^2 values.
r.squaredGLMM(CTM.model.m1)## R2m R2c
## [1,] 0.4683516 0.508793# Surprisingly reasonable. 47% of response explained by offspring temperature acclimation
# Calculating p-values
CTM.dwdata1 <- data.frame(term = c(names(summary(CTM.model.m1)$coefficients$fixed)), estimate = c(summary(CTM.model.m1)$coefficients$fixed[1:2]), std.error = c(summary(CTM.model.m1)$tTable[, 2]), statistic = c(summary(CTM.model.m1)$tTable[, 4]))
p.values <- c()
for (i in 1:nrow(CTM.dwdata1)) {
p.values[i] <- pt(abs(CTM.dwdata1$statistic[i]), df = 8, lower.tail = FALSE)
}
CTM.dwdata1$p.values <- round(p.values, 3)
CTM.dwdata1$term <- c("(Intercept)", "Offspring Acclimation Temperature")
# Plotting fitted model
CTM.Fits <- ggplot(
data = data.frame(
Pred = c(predict(CTM.model.m1)), Resp = c(CTM.NMS.Data$CTM),
Model = "Mod1"
),
aes(x = Pred, y = Resp, fill = Model, colour = Model)
) +
theme_bw() +
my.theme +
geom_point(pch = 21, size = 3, alpha = 0.5) +
geom_smooth(method = "lm") +
xlab("Predicted Critical Thermal Maximum") +
ylab("Measured Critical \n Thermal Maximum") +
scale_fill_manual(values = c("cornflowerblue")) +
scale_colour_manual(values = c("black")) +
theme(legend.position = "none") +
xlim(c(28.4, 29.2))
print(CTM.Fits)
ggsave("Critical Thermal Maximum Fit Comparisons.jpeg", dpi = 1000, limitsize = TRUE)
CTM.Models <- CTM.dwdata1 %>%
mutate(model = "Model 1")Testing model assumptions for resting metabic rate: model 1.
RMR.m1 <- lmer(log(RMR) ~ log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID), REML = FALSE, data = RMR.NMS.Data)
# Checking for collineary between predictors.
car::vif(RMR.m1) # very low evidence of collinearity between predictors.## log(Mass) Mother.Acc
## 1.011368 1.011368# R2?
r.squaredGLMM(RMR.m1) # low marginal R2 (0.134)## R2m R2c
## [1,] 0.1339877 0.1873523require(r2glmm)
r2beta(model = RMR.m1, partial = TRUE, method = "nsj") # Mass partial R2 > that of maternal acclimation temperatures## Effect Rsq upper.CL lower.CL
## 1 Model 0.134 0.235 0.061
## 2 log(Mass) 0.095 0.185 0.031
## 3 Mother.Acc2 0.060 0.140 0.012summary(RMR.m1)$varcor## Groups Name Std.Dev.
## Mother.ID (Intercept) 0.00000
## Father.ID (Intercept) 0.11360
## Residual 0.44331# Minimal variance explained by random effects (which explains singular fit), but retained for comparability. Diagnosing residuals
{
Res.Alone <- ggplot(
data = data.frame("Res" = residuals(RMR.m1, type = "deviance")),
aes(x = 1:nrow(RMR.NMS.Data), y = Res)
) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) +
xlab("Row Number") +
ylab("Deviance Residuals")
Res.by.Fit <- ggplot(
data = data.frame("Res" = residuals(RMR.m1), "Fit" = predict(RMR.m1)),
aes(x = Fit, y = Res)
) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) +
xlab("Fitted Values") +
ylab("Deviance Residuals")
Res.by.Resp <- ggplot(data = data.frame(
"Res" = residuals(RMR.m1), "Resp" =
log(RMR.NMS.Data$RMR)
), aes(x = Resp, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) +
xlab("Resting Metabolic Rate") +
ylab("Deviance Residuals")
QQNorm <- ggplot(
data.frame("Res" = residuals(RMR.m1), "Fit" = predict(RMR.m1)),
aes(sample = Res)
) +
stat_qq(colour = "gold") +
stat_qq_line() +
my.theme +
theme_bw()
}
grid.arrange(Res.Alone, Res.by.Fit, Res.by.Resp, QQNorm, nrow = 2, top = "Resting Metabolic Rate Residuals")
## Residuals posisbly non-normal; slight high-end taper. Checking confidence intervals
qqPlot(residuals(RMR.m1, type = "deviance"), ylab = "Deviance Residuals")
## 83 81
## 68 66# Taper is within confidence intervals, however. Using a histogram to confirm this speculation below.
Resid.hist <- ggplot(
data.frame(Res = c(residuals(RMR.m1, type = "deviance")), Model = "Model1"),
aes(Res, fill = Model, colour = Model)
) +
geom_histogram(alpha = 0.5, bins = 20) +
theme_bw() +
my.theme +
xlab("Deviance Residuals") +
ylab("Count") +
scale_fill_manual(
values =
c("cornflowerblue")
) +
scale_colour_manual(values = c("black")) +
theme(legend.position = "none") +
geom_density(alpha = 0.5, aes(y = 1.0 / 9 * ..count..))
print(Resid.hist)
# Yes, reasonably normal distribution. Plotting residuals by predictors.
{
Res.by.Mass <- ggplot(data = data.frame(
"Mass" = RMR.NMS.Data$Mass,
"Res" = residuals(RMR.m1)
), aes(x = Mass, y = Res, colour = as.factor(RMR.NMS.Data$Mother.Acc))) +
theme_bw() +
geom_point(size = 3, alpha = 0.5) +
my.theme +
xlab("Mass (g)") +
ylab("Deviance Residuals") +
scale_colour_manual(values = c("royalblue", "gold2")) +
guides(colour = guide_legend(title = "Dam Acclimation \nTemperature"))
Res.by.Dam <- ggplot(data = data.frame(
"Dam" = as.factor(RMR.NMS.Data$Mother.Acc),
"Res" = residuals(RMR.m1)
), aes(x = Dam, y = Res)) +
theme_bw() +
geom_boxplot(fill = "royalblue") +
my.theme +
xlab("Dam Acclimation Temperature (Degrees C)") +
ylab("Deviance Residuals")
}
grid.arrange(Res.by.Mass, Res.by.Dam, nrow = 2, top = "Resting Metabolic Rate Residuals")
# Homoskedastic. Calculating p-values.
dwdata.RMR.NMS.m1 <- data.frame(
term = c(rownames(summary(RMR.m1)$coefficients)),
estimate = c(summary(RMR.m1)$coefficients[1:3]),
std.error = c(summary(RMR.m1)$coefficients[4:6]),
statistic = c(summary(RMR.m1)$coefficients[7:9])
)
p.values.m1 <- c()
for (i in 1:nrow(dwdata.RMR.NMS.m1)) {
p.values.m1[i] <- pt(abs(dwdata.RMR.NMS.m1$statistic[i]), df = 8, lower.tail = FALSE)
}
dwdata.RMR.NMS.m1$p.values <- round(p.values.m1, 3)
dwdata.RMR.NMS.m1$term <- c("(Intercept)", "log Mass", "Dam Acclimation Temperature")Testing model assumptions for resting metabic rate: model 2.
RMR.m2 <- lmer(log(RMR) ~ log(Mass) + Offspring.Acc + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID),
REML = FALSE, data = RMR.NMS.Data
)
car::vif(RMR.m2) # again, low evidence of collinearity.## log(Mass) Offspring.Acc Mother.Acc
## 1.170197 1.158577 1.018361r.squaredGLMM(RMR.m2)## R2m R2c
## [1,] 0.141776 0.2053687# Checking partial-R2 values with Nakagawa and Schielzeth (2013) method for estimating DF.
r2beta(model = RMR.m2, partial = TRUE, method = "nsj") # Mass R^2 has greatest contribution.## Effect Rsq upper.CL lower.CL
## 1 Model 0.142 0.248 0.071
## 2 log(Mass) 0.101 0.193 0.035
## 4 Mother.Acc2 0.063 0.144 0.013
## 3 Offspring.Acc2 0.005 0.045 0.000# Checking variance explained by randon intercepts
summary(RMR.m2)$varcor## Groups Name Std.Dev.
## Mother.ID (Intercept) 0.00000
## Father.ID (Intercept) 0.12490
## Residual 0.44152# Again, maternal ID alone explaining zero variance. Plotting residuals by fitteds
{
Res.Alone <- ggplot(
data = data.frame("Res" = residuals(RMR.m2, type = "deviance")),
aes(x = 1:nrow(RMR.NMS.Data), y = Res)
) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) +
xlab("Row Number") +
ylab("Deviance Residuals")
Res.by.Fit <- ggplot(data = data.frame(
"Res" = residuals(RMR.m2, type = "deviance"),
"Fit" = predict(RMR.m2)
), aes(x = Fit, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) +
xlab("Fitted Values") +
ylab("Deviance Residuals")
Res.by.Resp <- ggplot(data = data.frame(
"Res" = residuals(RMR.m2, type = "deviance"),
"Resp" = log(RMR.NMS.Data$RMR)
), aes(x = Resp, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) +
xlab("Resting Metabolic Rate") +
ylab("Deviance Residuals")
QQNorm <- ggplot(data.frame(
"Res" = residuals(RMR.m2, type = "deviance"),
"Fit" = predict(RMR.m2)
), aes(sample = Res)) +
stat_qq(colour = "gold") +
stat_qq_line() +
my.theme +
theme_bw()
}
grid.arrange(Res.Alone, Res.by.Fit, Res.by.Resp, QQNorm, nrow = 2, top = "Resting Metabolic Rate Residuals")
## Residuals similar to model 1.
qqPlot(residuals(RMR.m2, type = "deviance"), ylab = "Deviance Residuals")
## 83 81
## 68 66# Again, generally within confidence intervals.
{
Res.by.Mass <- ggplot(data = data.frame(
"Mass" = RMR.NMS.Data$Mass,
"Res" = residuals(RMR.m2)
), aes(x = Mass, y = Res, colour = as.factor(RMR.NMS.Data$Mother.Acc))) +
theme_bw() +
geom_point(size = 3, alpha = 0.5) +
my.theme +
xlab("Mass (g)") +
ylab("Deviance Residuals") +
scale_colour_manual(values = c("royalblue", "gold2")) +
guides(colour = guide_legend(title = "Dam Acclimation \nTemperature"))
Res.by.Dam <- ggplot(data = data.frame(
"Dam" = as.factor(RMR.NMS.Data$Mother.Acc),
"Res" = residuals(RMR.m2)
), aes(x = Dam, y = Res)) +
theme_bw() +
geom_boxplot(fill = "royalblue") +
my.theme.diag +
xlab("Dam Acclimation Temperature \n (Degrees C)") +
ylab("Deviance Residuals")
Res.by.Off <- ggplot(data = data.frame(
"Offspring" = as.factor(RMR.NMS.Data$Offspring.Acc),
"Res" = residuals(RMR.m2)
), aes(x = Offspring, y = Res)) +
theme_bw() +
geom_boxplot(fill = "gold2") +
my.theme.diag +
xlab("Dam Acclimation Temperature \n (Degrees C)") +
ylab("Deviance Residuals")
}
layout <- rbind(c(1, 1), c(2, 3))
grid.arrange(Res.by.Mass, Res.by.Dam, Res.by.Off, nrow = 2, top = "Resting Metabolic Rate Residuals", layout_matrix = layout)
# Again, homoskedastic. Calculating p-values.
dwdata.RMR.NMS.m2 <- data.frame(
term = c(rownames(summary(RMR.m2)$coefficients)),
estimate = c(summary(RMR.m2)$coefficients[1:4]),
std.error = c(summary(RMR.m2)$coefficients[5:8]),
statistic = c(summary(RMR.m2)$coefficients[9:12])
)
p.values.m2 <- c()
for (i in 1:nrow(dwdata.RMR.NMS.m2)) {
p.values.m2[i] <- pt(abs(dwdata.RMR.NMS.m2$statistic[i]), df = 8, lower.tail = FALSE)
}
dwdata.RMR.NMS.m2$p.values <- round(p.values.m2, 3)
dwdata.RMR.NMS.m2$term <- c("(Intercept)", "log Mass", "Offspring Acclimation Temperature", "Dam Acclimation Temperature")Testing model assumptions for resting metabic rate: model 3.
RMR.m3 <- lmer(log(RMR) ~ log(Mass) + Offspring.Acc + Mother.Acc + Offspring.Acc:Mother.Acc +
(1 | Mother.ID) + (1 | Father.ID), REML = FALSE, data = RMR.NMS.Data)
summary(RMR.m3)$varcor## Groups Name Std.Dev.
## Mother.ID (Intercept) 0.0083881
## Father.ID (Intercept) 0.1284166
## Residual 0.4390322# As before, maternal ID alone explains nearly zero variance, but estimate is possible (thus, no singular fit error). Plotting residuals.
{
Res.Alone <- ggplot(
data = data.frame("Res" = residuals(RMR.m3, type = "deviance")),
aes(x = 1:nrow(RMR.NMS.Data), y = Res)
) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) +
xlab("Row Number") +
ylab("Deviance Residuals")
Res.by.Fit <- ggplot(data = data.frame(
"Res" = residuals(RMR.m3, type = "deviance"),
"Fit" = predict(RMR.m3)
), aes(x = Fit, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) +
xlab("Fitted Values") +
ylab("Deviance Residuals")
Res.by.Resp <- ggplot(data = data.frame(
"Res" = residuals(RMR.m3, type = "deviance"),
"Resp" = log(RMR.NMS.Data$RMR)
), aes(x = Resp, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) +
xlab("Resting Metabolic Rate") +
ylab("Deviance Residuals")
QQNorm <- ggplot(data.frame(
"Res" = residuals(RMR.m3, type = "deviance"),
"Fit" = predict(RMR.m3)
), aes(sample = Res)) +
stat_qq(colour = "gold") +
stat_qq_line() +
my.theme +
theme_bw()
}
grid.arrange(Res.Alone, Res.by.Fit, Res.by.Resp, QQNorm, nrow = 2, top = "Resting Metabolic Rate Residuals")
## Very subtle tails in residuals...
qqPlot(residuals(RMR.m3, type = "deviance"), ylab = "Deviance Residuals")
## 83 81
## 68 66## But largely within confidence intervals.
Resid.hist <- ggplot(
data.frame(Res = c(residuals(RMR.m3, type = "deviance")), Model = "Model1"),
aes(Res, fill = Model, colour = Model)
) +
geom_histogram(alpha = 0.5, bins = 20) +
theme_bw() +
my.theme +
xlab("Deviance Residuals") +
ylab("Count") +
scale_fill_manual(
values =
c("firebrick")
) +
scale_colour_manual(values = c("black")) +
theme(legend.position = "none") +
geom_density(alpha = 0.5, aes(y = 0.5 / 5 * ..count..))
print(Resid.hist)
# Very slightly bimodal, but acceptable. Assessing residuals across predictors.
{
Res.by.Mass <- ggplot(data = data.frame(
"Mass" = RMR.NMS.Data$Mass,
"Res" = residuals(RMR.m2)
), aes(x = Mass, y = Res, colour = as.factor(RMR.NMS.Data$Mother.Acc))) +
theme_bw() +
geom_point(size = 3, alpha = 0.5) +
my.theme.diag +
xlab("Mass (g)") +
ylab("Deviance Residuals") +
scale_colour_manual(values = c("royalblue", "gold2")) +
guides(colour = guide_legend(title = "Dam Acclimation Temperature")) +
theme(legend.position = "top")
Res.by.Dam <- ggplot(data = data.frame(
"Dam" = as.factor(RMR.NMS.Data$Mother.Acc),
"Res" = residuals(RMR.m2)
), aes(x = Dam, y = Res)) +
theme_bw() +
geom_boxplot(fill = "royalblue") +
my.theme.diag +
xlab("Dam Acclimation Temperature \n (Degrees C)") +
ylab("Deviance Residuals")
Res.by.Off <- ggplot(data = data.frame(
"Offspring" = as.factor(RMR.NMS.Data$Offspring.Acc),
"Res" = residuals(RMR.m2)
), aes(x = Offspring, y = Res)) +
theme_bw() +
geom_boxplot(fill = "gold2") +
my.theme.diag +
xlab("Dam Acclimation Temperature \n (Degrees C)") +
ylab("Deviance Residuals")
Off.by.Dam <- ggplot(data = data.frame(
"Dam" = as.factor(RMR.NMS.Data$Mother.Acc),
"Res" = residuals(RMR.m2)
), aes(x = Dam, y = Res, fill = as.factor(RMR.NMS.Data$Offspring.Acc))) +
theme_bw() +
geom_boxplot() +
my.theme.diag +
scale_fill_manual(values = c("firebrick", "coral")) +
xlab("Dam Acclimation Temperature (Degrees C)") +
ylab("Deviance Residuals") +
guides(fill = guide_legend(title = "Offspring Acclimation \n Temperature")) +
theme(legend.position = "top")
}
layout <- rbind(c(1, 1, 2, 3), c(1, 1, 4, 4))
grid.arrange(Res.by.Mass, Res.by.Dam, Res.by.Off, Off.by.Dam, nrow = 2, top = "Resting Metabolic Rate Residuals", layout_matrix = layout)
# Appear acceptably homoskedastic. Calculating p-values.
dwdata.RMR.NMS.m3 <- data.frame(
term = c(rownames(summary(RMR.m3)$coefficients)),
estimate = c(summary(RMR.m3)$coefficients[1:5]),
std.error = c(summary(RMR.m3)$coefficients[6:10]),
statistic = c(summary(RMR.m3)$coefficients[11:15])
)
p.values.m3 <- c()
for (i in 1:nrow(dwdata.RMR.NMS.m3)) {
p.values.m3[i] <- pt(abs(dwdata.RMR.NMS.m3$statistic[i]), df = 8, lower.tail = FALSE)
}
dwdata.RMR.NMS.m3$p.values <- round(p.values.m3, 3)
dwdata.RMR.NMS.m3$term <- c("(Intercept)", "log Mass", "Offspring Acclimation Temperature", "Dam Acclimation Temperature", "Dam Acclimation Temperature :\n Offspring Acclimation Temperature")
## Comparing models by plot.
dwdata.RMR.NMS.m1 <- dwdata.RMR.NMS.m1 %>%
mutate(model = "Model 1")
dwdata.RMR.NMS.m2 <- dwdata.RMR.NMS.m2 %>%
mutate(model = "Model 2")
dwdata.RMR.NMS.m3 <- dwdata.RMR.NMS.m3 %>%
mutate(model = "Model 3")
RMR.NMS.Models <- rbind(dwdata.RMR.NMS.m1, dwdata.RMR.NMS.m2, dwdata.RMR.NMS.m3)
RMR.NMS.Models$estimate <- RMR.NMS.Models$estimate * 5
RMR.NMS.Models$std.error <- RMR.NMS.Models$std.error * 5
RMR.NMS.Models %>%
dwplot(style = "distribution", show_intercept = FALSE) + theme_bw() +
theme(legend.justification = c(-0.1, -0.1), legend.position = c(0.7, 0.01)) + my.theme.dw +
xlab("Coefficient") + ylab("") + geom_vline(xintercept = 0, colour = "black", linetype = 2) +
scale_fill_manual(values = c("mediumseagreen", "mediumorchid", "cornflowerblue")) +
scale_colour_manual(values = c("black", "black", "black")) + scale_x_continuous(
limits = c(-5.0, 5.0),
breaks = c(-5.0, -2.5, 0, 2.5, 5.0), label = c("-1.00", "-0.50", "0", "0.50", "1.00")
)
ggsave("log Resting Metabolic Rate Model Comparisons.jpeg", dpi = 1000, limitsize = TRUE)
RMR.NMS.Models$estimate <- RMR.NMS.Models$estimate / 5
RMR.NMS.Models$std.error <- RMR.NMS.Models$std.error / 5
RMR.NMS.Fit.m1 <- data.frame("Response" = c(log(RMR.NMS.Data$RMR)), "Fit" = c(predict(RMR.m1)))
RMR.NMS.Fit.m2 <- data.frame("Response" = c(log(RMR.NMS.Data$RMR)), "Fit" = c(predict(RMR.m2)))
RMR.NMS.Fit.m3 <- data.frame("Response" = c(log(RMR.NMS.Data$RMR)), "Fit" = c(predict(RMR.m3)))
RMR.NMS.Fit.m1 <- RMR.NMS.Fit.m1 %>%
mutate(model = "Model 1")
RMR.NMS.Fit.m2 <- RMR.NMS.Fit.m2 %>%
mutate(model = "Model 2")
RMR.NMS.Fit.m3 <- RMR.NMS.Fit.m3 %>%
mutate(model = "Model 3")
RMR.NMS.Fit.All <- rbind(RMR.NMS.Fit.m1, RMR.NMS.Fit.m2, RMR.NMS.Fit.m3)
RMR.NMS.Fits <- ggplot(data = RMR.NMS.Fit.All, aes(
x = Fit, y = Response, fill = model, colour = model,
linetype = model
)) +
theme_bw() +
my.theme +
geom_point(pch = 21, size = 3, alpha = 0.5) +
geom_smooth(method = "lm") +
xlab("Predicted log Resting Metabolic Rate \n (log mg O2/h)") +
ylab("Measured log Resting \n Metabolic Rate \n (log mg O2/h)") +
scale_fill_manual(values = c(
"mediumseagreen",
"mediumorchid", "cornflowerblue"
)) +
scale_colour_manual(values = c(
"black",
"black", "black"
)) +
theme(legend.title = element_blank())
print(RMR.NMS.Fits)
ggsave("log Resting Metabolic Rate Model Fit Comparisons.jpeg", dpi = 1000, limitsize = TRUE)Next, testing model assumptions for peak metabic rate: model 1.
MMR.m1 <- lmer(log(MMR) ~ log(Mass) + (1 | Mother.ID) + (1 | Father.ID), REML = FALSE, data = MMR.NMS.Data)
# R^2 values?
r.squaredGLMM(MMR.m1) # Reasonable variance explained## R2m R2c
## [1,] 0.3692743 0.3719902summary(MMR.m1)$varcor # Maternal effect non-detectable (as random intercept)## Groups Name Std.Dev.
## Mother.ID (Intercept) 0.000000
## Father.ID (Intercept) 0.016496
## Residual 0.250845{
Res.Alone <- ggplot(
data = data.frame("Res" = residuals(MMR.m1, type = "deviance")),
aes(x = 1:nrow(MMR.NMS.Data), y = Res)
) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) +
xlab("Row Number") +
ylab("Deviance Residuals")
Res.by.Fit <- ggplot(data = data.frame(
"Res" = residuals(MMR.m1, type = "deviance"),
"Fit" = predict(MMR.m1)
), aes(x = Fit, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) +
xlab("Fitted Values") +
ylab("Deviance Residuals")
Res.by.Resp <- ggplot(data = data.frame(
"Res" = residuals(MMR.m1, type = "deviance"),
"Resp" = MMR.NMS.Data$MMR
), aes(x = Resp, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) +
xlab("Resting Metabolic Rate") +
ylab("Deviance Residuals")
QQNorm <- ggplot(data.frame(
"Res" = residuals(MMR.m1, type = "deviance"),
"Fit" = predict(MMR.m1)
), aes(sample = Res)) +
stat_qq(colour = "gold") +
stat_qq_line() +
my.theme +
theme_bw()
}
grid.arrange(Res.Alone, Res.by.Fit, Res.by.Resp, QQNorm, nrow = 2, top = "Peak Metabolic Rate Residuals")
# Wow! One very destinct outlier. Checking individual.
log(MMR.NMS.Data[which(residuals(MMR.m1, type = "deviance") < -2.0), "MMR"])## [1] -2.149864MMR.NMS.Data %>%
summarise(Mean_MMR = mean(log(MMR), na.rm = T), "LC" = mean(log(MMR), na.rm = T) - 4 * sd(log(MMR), na.rm = T))## Mean_MMR LC
## 1 0.1334902 -1.13047# Remarkably low peak metabolic rate and absolute aerobic scope. Removing from both datasets.
MMR.NMS.Data <- MMR.NMS.Data[-c(which(residuals(MMR.m1, type = "deviance") < -2.0)), ]Recalculting AICc values for peak metabolic rate.
MMR.model <- lmer(log(MMR) ~ Mother.Acc * Father.Acc * Offspring.Acc + log(Mass) + (1 | Mother.ID) + (1 | Father.ID), REML = FALSE, data = MMR.NMS.Data)
AIC.Results.MMR <- subset(dredge(MMR.model, rank = "AICc", extra = list(R2 = function(x) r.squaredGLMM(x))), delta < 2)
Select.Models.MMR <- c()
Delta.MMR <- c()
Response <- c()
Marg.R2 <- c()
Cond.R2 <- c()
for (i in 1:length(attr(AIC.Results.MMR, "model.calls"))) {
Select.Models.MMR[i] <- gsub(".*~ ", "", paste(attr(AIC.Results.MMR, "model.calls")[[i]])[2])
Delta.MMR[i] <- AIC.Results.MMR$delta[i]
Response[i] <- gsub(" ~.*", "", paste(attr(AIC.Results.MMR, "model.calls")[[i]])[2])
Marg.R2[i] <- AIC.Results.MMR$R21[i]
Cond.R2[i] <- AIC.Results.MMR$R22[i]
if (i == length(attr(AIC.Results.MMR, "model.calls"))) {
Table <- data.frame(
"Reponse.Type" = "Mass Independant", "Response" = Response,
"Delta.AIC" = Delta.MMR, "Marginal.R2" = Marg.R2, "Conditional.R2" = Cond.R2,
"Equation" = Select.Models.MMR
)
assign(paste("AIC.Selection", "MMR", sep = "_"), Table)
}
}
AIC.Selection.NMS <- rbind(AIC.Selection_NMS_1, AIC.Selection_MMR, AIC.Selection_NMS_3)
View(AIC.Selection.NMS)
write.csv(AIC.Selection.NMS, "Model_AIC_Final.csv", row.names = FALSE)Reassessing model assumptions for peak metabic rate: model 1.
MMR.m1 <- lmer(log(MMR) ~ log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID), REML = FALSE, data = MMR.NMS.Data)
# Checking R^2 values
r.squaredGLMM(MMR.m1) # Considerable improvement## R2m R2c
## [1,] 0.5107389 0.5241443r2beta(model = MMR.m1, partial = TRUE, method = "nsj") # Mass explains 50.5% of variation.## Effect Rsq upper.CL lower.CL
## 1 Model 0.511 0.589 0.431
## 2 log(Mass) 0.505 0.583 0.424
## 3 Offspring.Acc2 0.031 0.090 0.002summary(MMR.m1)$AICtab[1]## AIC
## -82.4201# AIC dramatically improved. Checking contribution of random intercepts
summary(MMR.m1)$varcor## Groups Name Std.Dev.
## Mother.ID (Intercept) 0.000000
## Father.ID (Intercept) 0.032587
## Residual 0.194153# Maternal ID again explaining no residual variance. Residual diagnostics:
{
Res.Alone <- ggplot(
data = data.frame("Res" = residuals(MMR.m1, type = "deviance")),
aes(x = 1:nrow(MMR.NMS.Data), y = Res)
) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) +
xlab("Row Number") +
ylab("Deviance Residuals")
Res.by.Fit <- ggplot(data = data.frame(
"Res" = residuals(MMR.m1, type = "deviance"),
"Fit" = predict(MMR.m1)
), aes(x = Fit, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) +
xlab("Fitted Values") +
ylab("Deviance Residuals")
Res.by.Resp <- ggplot(data = data.frame(
"Res" = residuals(MMR.m1, type = "deviance"),
"Resp" = MMR.NMS.Data$MMR
), aes(x = Resp, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) +
xlab("Resting Metabolic Rate") +
ylab("Deviance Residuals")
QQNorm <- ggplot(data.frame(
"Res" = residuals(MMR.m1, type = "deviance"),
"Fit" = predict(MMR.m1)
), aes(sample = Res)) +
stat_qq(colour = "gold") +
stat_qq_line() +
my.theme +
theme_bw()
}
grid.arrange(Res.Alone, Res.by.Fit, Res.by.Resp, QQNorm, nrow = 2, top = "Peak Metabolic Rate Residuals")
qqPlot(residuals(MMR.m1, type = "deviance"), ylab = "Deviance Residuals")
## 198 47
## 194 47# Residuals subtly non-normal. Double-checking with a histogram
Resid.hist <- ggplot(
data.frame(Res = c(residuals(MMR.m1, type = "deviance")), Model = "Model1"),
aes(Res, fill = Model, colour = Model)
) +
geom_histogram(alpha = 0.5, bins = 20) +
theme_bw() +
my.theme +
xlab("Deviance Residuals") +
ylab("Count") +
scale_fill_manual(
values =
c("black")
) +
scale_colour_manual(values = c("navajowhite")) +
theme(legend.position = "none") +
geom_density(alpha = 0.5, aes(y = 0.2 / 5 * ..count..))
print(Resid.hist)
# Not bad, although appearing slightly bimodal. This could, simply, be a consequence of limited data. Checking for heteroskedasticity across predictor variables nonetheless.
{
Res.by.Mass <- ggplot(data = data.frame(
"Mass" = RMR.NMS.Data$Mass,
"Res" = residuals(RMR.m2)
), aes(x = Mass, y = Res, colour = as.factor(RMR.NMS.Data$Offspring.Acc))) +
theme_bw() +
geom_point(size = 3, alpha = 0.5) +
my.theme.diag +
xlab("Mass (g)") +
ylab("Deviance Residuals") +
scale_colour_manual(values = c("royalblue", "gold2")) +
guides(colour = guide_legend(title = "Offspring Acclimation Temperature")) +
theme(legend.position = "top")
Res.by.Off <- ggplot(data = data.frame(
"Off" = as.factor(RMR.NMS.Data$Offspring.Acc),
"Res" = residuals(RMR.m2)
), aes(x = Off, y = Res)) +
theme_bw() +
geom_boxplot(fill = "royalblue") +
my.theme.diag +
xlab("Offspring Acclimation Temperature (Degrees C)") +
ylab("Deviance Residuals")
QQNorm.Diag <- ggplot(
data.frame(
"Res" = residuals(MMR.m1, type = "deviance"),
"Fit" = predict(MMR.m1), "Off" = as.factor(MMR.NMS.Data$Offspring.Acc)
),
aes(sample = Res, colour = Off)
) +
stat_qq() +
stat_qq_line() +
my.theme +
theme_bw() +
scale_colour_manual(values = c("black", "thistle")) +
theme(legend.position = "top") +
guides(colour = guide_legend(title = "Offspring Acclimation Temperature"))
}
layout <- rbind(c(1, 1, 2, 2), c(NA, 3, 3, NA))
grid.arrange(Res.by.Mass, Res.by.Off, QQNorm.Diag, nrow = 2, top = "Peak Metabolic Rate Residuals", layout_matrix = layout)
# No clear heteroskedasticity, nor difference in sample size. Assessing with a histogram.
Resid.hist <- ggplot(
data.frame(
Res = c(residuals(MMR.m1, type = "deviance")),
"Offspring" = as.factor(MMR.NMS.Data$Offspring.Acc)
),
aes(Res, fill = Offspring, colour = Offspring)
) +
geom_histogram(alpha = 0.5, bins = 20) +
theme_bw() +
my.theme +
xlab("Deviance Residuals") +
ylab("Count") +
scale_fill_manual(
values =
c("black", "mediumorchid")
) +
scale_colour_manual(values = c("white", "black")) +
theme(legend.position = "none") +
geom_density(alpha = 0.5, aes(y = 0.2 / 4 * ..count..))
print(Resid.hist)
## Equivalent data distributions. Continuing on to produce p-values.
dwdata.MMR.NMS.m1 <- data.frame(
term = c(rownames(summary(MMR.m1)$coefficients)),
estimate = c(summary(MMR.m1)$coefficients[1:3]),
std.error = c(summary(MMR.m1)$coefficients[4:6]),
statistic = c(summary(MMR.m1)$coefficients[7:9])
)
p.values.m1 <- c()
for (i in 1:nrow(dwdata.MMR.NMS.m1)) {
p.values.m1[i] <- pt(abs(dwdata.MMR.NMS.m1$statistic[i]), df = 8, lower.tail = FALSE)
}
dwdata.MMR.NMS.m1$p.values <- round(p.values.m1, 3)
dwdata.MMR.NMS.m1$term <- c("(Intercept)", "log Mass", "Offspring Acclimation Temperature")Testing model assumptions for peak metabic rate: model 2.
MMR.m2 <- lmer(log(MMR) ~ log(Mass) + Offspring.Acc + Father.Acc + (1 | Mother.ID) + (1 | Father.ID), REML = FALSE, data = MMR.NMS.Data)
# Contribution of random intercepts?
summary(MMR.m2)$varcor## Groups Name Std.Dev.
## Mother.ID (Intercept) 0.000000
## Father.ID (Intercept) 0.027424
## Residual 0.194068# Maternal and paternal ID explain very little residual variance.
{
Res.Alone <- ggplot(
data = data.frame("Res" = residuals(MMR.m2, type = "deviance")),
aes(x = 1:nrow(MMR.NMS.Data), y = Res)
) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) +
xlab("Row Number") +
ylab("Deviance Residuals")
Res.by.Fit <- ggplot(data = data.frame(
"Res" = residuals(MMR.m2, type = "deviance"),
"Fit" = predict(MMR.m2)
), aes(x = Fit, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) +
xlab("Fitted Values") +
ylab("Deviance Residuals")
Res.by.Resp <- ggplot(data = data.frame(
"Res" = residuals(MMR.m2, type = "deviance"),
"Resp" = MMR.NMS.Data$MMR
), aes(x = Resp, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) +
xlab("Resting Metabolic Rate") +
ylab("Deviance Residuals")
QQNorm <- ggplot(data.frame(
"Res" = residuals(MMR.m2, type = "deviance"),
"Fit" = predict(MMR.m2)
), aes(sample = Res)) +
stat_qq(colour = "gold") +
stat_qq_line() +
my.theme +
theme_bw()
}
grid.arrange(Res.Alone, Res.by.Fit, Res.by.Resp, QQNorm, nrow = 2, top = "Peak Metabolic Rate Residuals")
# Residuals appear fair, although straying from normality as observed before.
qqPlot(residuals(MMR.m2, type = "deviance"), ylab = "Deviance Residuals")
## 198 47
## 194 47Resid.hist <- ggplot(
data.frame(Res = c(residuals(MMR.m2, type = "deviance")), Model = "Model1"),
aes(Res, fill = Model, colour = Model)
) +
geom_histogram(alpha = 0.5, bins = 20) +
theme_bw() +
my.theme +
xlab("Deviance Residuals") +
ylab("Count") +
scale_fill_manual(
values =
c("cornflowerblue")
) +
scale_colour_manual(values = c("black")) +
theme(legend.position = "none") +
geom_density(alpha = 0.5, aes(y = 0.2 / 5 * ..count..))
print(Resid.hist)
# Limited concern in histogram. Assessing parameter-specific variance.
{
Res.by.Mass <- ggplot(data = data.frame(
"Mass" = MMR.NMS.Data$Mass,
"Res" = residuals(MMR.m2, type = "deviance")
), aes(
x = Mass, y = Res,
colour = as.factor(MMR.NMS.Data$Offspring.Acc)
)) +
theme_bw() +
geom_point(size = 3, alpha = 0.5) +
my.theme.diag +
xlab("Mass (g)") +
ylab("Deviance Residuals") +
scale_colour_manual(values = c("royalblue", "gold2")) +
guides(colour = guide_legend(title = "Offspring Acclimation Temperature")) +
theme(legend.position = "top")
Res.by.Off <- ggplot(data = data.frame(
"Off" = as.factor(MMR.NMS.Data$Offspring.Acc),
"Res" = residuals(MMR.m2, type = "deviance")
), aes(x = Off, y = Res)) +
theme_bw() +
geom_boxplot(fill = "royalblue") +
my.theme.diag +
xlab("Offspring Acclimation Temperature (Degrees C)") +
ylab("Deviance Residuals")
Res.by.Pat <- ggplot(data = data.frame(
"Pat" = as.factor(MMR.NMS.Data$Father.Acc),
"Res" = residuals(MMR.m2, type = "deviance")
), aes(x = Pat, y = Res)) +
theme_bw() +
geom_boxplot(fill = "mediumseagreen") +
my.theme.diag +
xlab("Sire Acclimation Temperature (Degrees C)") +
ylab("Deviance Residuals")
QQNorm.by.Off <- ggplot(
data.frame(
"Res" = residuals(MMR.m2, type = "deviance"),
"Fit" = predict(MMR.m2), "Off" = as.factor(MMR.NMS.Data$Offspring.Acc)
),
aes(sample = Res, colour = Off)
) +
stat_qq() +
stat_qq_line() +
my.theme +
theme_bw() +
scale_colour_manual(values = c("mediumseagreen", "mediumorchid")) +
theme(legend.position = "top") +
guides(colour = guide_legend(title = "Offspring Acclimation Temperature"))
QQNorm.by.Pat <- ggplot(
data.frame(
"Res" = residuals(MMR.m2, type = "deviance"),
"Fit" = predict(MMR.m2), "Pat" = as.factor(MMR.NMS.Data$Father.Acc)
),
aes(sample = Res, colour = Pat)
) +
stat_qq() +
stat_qq_line() +
my.theme +
theme_bw() +
scale_colour_manual(values = c("royalblue", "gold2")) +
theme(legend.position = "top") +
guides(colour = guide_legend(title = "Sire Acclimation Temperature"))
}
grid.arrange(Res.by.Mass, Res.by.Off, Res.by.Pat, QQNorm.by.Off, QQNorm.by.Pat, nrow = 3, top = "Peak Metabolic Rate Residuals")
# Perhaps a slight heterogeneity across offspring acclimation temperature, but very minimal. Little reason for concern.
# Calculating p-values.
dwdata.MMR.NMS.m2 <- data.frame(
term = c(rownames(summary(MMR.m2)$coefficients)),
estimate = c(summary(MMR.m2)$coefficients[1:4]),
std.error = c(summary(MMR.m2)$coefficients[5:8]),
statistic = c(summary(MMR.m2)$coefficients[9:12])
)
p.values.m2 <- c()
for (i in 1:nrow(dwdata.MMR.NMS.m2)) {
p.values.m2[i] <- pt(abs(dwdata.MMR.NMS.m2$statistic[i]), df = 8, lower.tail = FALSE)
}
dwdata.MMR.NMS.m2$p.values <- round(p.values.m2, 3)
dwdata.MMR.NMS.m2$term <- c("(Intercept)", "log Mass", "Offspring Acclimation Temperature", "Sire Acclimation Temperature")Testing model assumptions for peak metabic rate: model 3.
MMR.m3 <- lmer(log(MMR) ~ log(Mass) + Offspring.Acc + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID), REML = FALSE, data = MMR.NMS.Data)
summary(MMR.m3)$varcor## Groups Name Std.Dev.
## Mother.ID (Intercept) 1.3465e-09
## Father.ID (Intercept) 3.3598e-02
## Residual 1.9388e-01# Paternal ID alone explaining residual variance. This pattern appears rather consistent.
{
Res.Alone <- ggplot(
data = data.frame("Res" = residuals(MMR.m3, type = "deviance")),
aes(x = 1:nrow(MMR.NMS.Data), y = Res)
) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) +
xlab("Row Number") +
ylab("Deviance Residuals")
Res.by.Fit <- ggplot(
data = data.frame("Res" = residuals(MMR.m3), "Fit" = predict(MMR.m3)),
aes(x = Fit, y = Res)
) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) +
xlab("Fitted Values") +
ylab("Deviance Residuals")
Res.by.Resp <- ggplot(data = data.frame(
"Res" = residuals(MMR.m3), "Resp" =
MMR.NMS.Data$MMR
), aes(x = Resp, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) +
xlab("Resting Metabolic Rate") +
ylab("Deviance Residuals")
QQNorm <- ggplot(
data.frame("Res" = residuals(MMR.m3), "Fit" = predict(MMR.m2)),
aes(sample = Res)
) +
stat_qq(colour = "gold") +
stat_qq_line() +
my.theme +
theme_bw()
}
grid.arrange(Res.Alone, Res.by.Fit, Res.by.Resp, QQNorm, nrow = 2, top = "Peak Metabolic Rate Residuals")
# Similar patterns to previous two models.
Resid.hist <- ggplot(
data.frame(Res = c(residuals(MMR.m3, type = "deviance")), Model = "Model1"),
aes(Res, fill = Model, colour = Model)
) +
geom_histogram(alpha = 0.5, bins = 20) +
theme_bw() +
my.theme +
xlab("Deviance Residuals") +
ylab("Count") +
scale_fill_manual(
values =
c("cornflowerblue")
) +
scale_colour_manual(values = c("black")) +
theme(legend.position = "none") +
geom_density(alpha = 0.5, aes(y = 0.2 / 5 * ..count..))
print(Resid.hist) # Little reason for concern. Plotting residuals by predictors.
{
Res.by.Mass <- ggplot(data = data.frame(
"Mass" = MMR.NMS.Data$Mass,
"Res" = residuals(MMR.m3, type = "deviance")
), aes(
x = Mass, y = Res,
colour = as.factor(MMR.NMS.Data$Offspring.Acc)
)) +
theme_bw() +
geom_point(size = 3, alpha = 0.5) +
my.theme.diag +
xlab("Mass (g)") +
ylab("Deviance Residuals") +
scale_colour_manual(values = c("royalblue", "gold2")) +
guides(colour = guide_legend(title = "Offspring Acclimation Temperature")) +
theme(legend.position = "top")
Res.by.Off <- ggplot(data = data.frame(
"Off" = as.factor(MMR.NMS.Data$Offspring.Acc),
"Res" = residuals(MMR.m3, type = "deviance")
), aes(x = Off, y = Res)) +
theme_bw() +
geom_boxplot(fill = "royalblue") +
my.theme.diag +
xlab("Offspring Acclimation Temperature (Degrees C)") +
ylab("Deviance Residuals")
Res.by.Mat <- ggplot(data = data.frame(
"Mat" = as.factor(MMR.NMS.Data$Mother.Acc),
"Res" = residuals(MMR.m3, type = "deviance")
), aes(x = Mat, y = Res)) +
theme_bw() +
geom_boxplot(fill = "mediumseagreen") +
my.theme.diag +
xlab("Dam Acclimation Temperature (Degrees C)") +
ylab("Deviance Residuals")
QQNorm.by.Off <- ggplot(
data.frame(
"Res" = residuals(MMR.m3, type = "deviance"),
"Fit" = predict(MMR.m3), "Off" = as.factor(MMR.NMS.Data$Offspring.Acc)
),
aes(sample = Res, colour = Off)
) +
stat_qq() +
stat_qq_line() +
my.theme +
theme_bw() +
scale_colour_manual(values = c("mediumseagreen", "mediumorchid")) +
theme(legend.position = "top") +
guides(colour = guide_legend(title = "Offspring Acclimation Temperature"))
QQNorm.by.Mat <- ggplot(
data.frame(
"Res" = residuals(MMR.m3, type = "deviance"),
"Fit" = predict(MMR.m3), "Mat" = as.factor(MMR.NMS.Data$Mother.Acc)
),
aes(sample = Res, colour = Mat)
) +
stat_qq() +
stat_qq_line() +
my.theme +
theme_bw() +
scale_colour_manual(values = c("royalblue", "gold2")) +
theme(legend.position = "top") +
guides(colour = guide_legend(title = "Dam Acclimation Temperature"))
}
grid.arrange(Res.by.Mass, Res.by.Off, Res.by.Mat, QQNorm.by.Off, QQNorm.by.Pat, nrow = 3, top = "Peak Metabolic Rate Residuals")
## Similar variance as above two models. Perhaps an interaction between mass and temperature categories
## could be considered? Calculating p-values.
dwdata.MMR.NMS.m3 <- data.frame(
term = c(rownames(summary(MMR.m3)$coefficients)),
estimate = c(summary(MMR.m3)$coefficients[1:4]),
std.error = c(summary(MMR.m3)$coefficients[5:8]),
statistic = c(summary(MMR.m3)$coefficients[9:12])
)
p.values.m3 <- c()
for (i in 1:nrow(dwdata.MMR.NMS.m3)) {
p.values.m3[i] <- pt(abs(dwdata.MMR.NMS.m3$statistic[i]), df = 8, lower.tail = FALSE)
}
dwdata.MMR.NMS.m3$p.values <- round(p.values.m3, 3)
dwdata.MMR.NMS.m3$term <- c("(Intercept)", "log Mass", "Offspring Acclimation Temperature", "Dam Acclimation Temperature")
## Comparing model outputs by plots
dwdata.MMR.NMS.m1 <- dwdata.MMR.NMS.m1 %>%
mutate(model = "Model 1")
dwdata.MMR.NMS.m2 <- dwdata.MMR.NMS.m2 %>%
mutate(model = "Model 2")
dwdata.MMR.NMS.m3 <- dwdata.MMR.NMS.m3 %>%
mutate(model = "Model 3")
MMR.NMS.Models <- rbind(dwdata.MMR.NMS.m1, dwdata.MMR.NMS.m2, dwdata.MMR.NMS.m3)
MMR.NMS.Models$estimate <- MMR.NMS.Models$estimate * 8
MMR.NMS.Models$std.error <- MMR.NMS.Models$std.error * 8
MMR.NMS.Models %>%
dwplot(style = "distribution", show_intercept = FALSE) + theme_bw() +
theme(
legend.justification = c(-0.1, -0.1),
legend.position = c(0, 0.01)
) + my.theme.dw + xlab("Coefficient") + ylab("") +
geom_vline(xintercept = 0, colour = "black", linetype = 2) +
scale_fill_manual(values = c("mediumseagreen", "mediumorchid", "cornflowerblue")) +
scale_colour_manual(values = c("black", "black", "black")) + scale_x_continuous(
limits = c(-8.0, 8.0),
breaks = c(-8.0, -4.0, 0, 4.0, 8.0), label = c("-1.00", "-0.50", "0", "0.50", "1.00")
)
ggsave("log Peak Metabolic Rate Model Comparisons.jpeg", dpi = 1000, limitsize = TRUE)
MMR.NMS.Models$estimate <- MMR.NMS.Models$estimate / 8
MMR.NMS.Models$std.error <- MMR.NMS.Models$std.error / 8
{
MMR.NMS.Fit.m1 <- data.frame(
"Response" = c(log(MMR.NMS.Data$MMR)), "Fit" = c(predict(MMR.m1)),
"Mass" = c(log(MMR.NMS.Data$Mass))
)
MMR.NMS.Fit.m2 <- data.frame(
"Response" = c(log(MMR.NMS.Data$MMR)), "Fit" = c(predict(MMR.m2)),
"Mass" = c(log(MMR.NMS.Data$Mass))
)
MMR.NMS.Fit.m3 <- data.frame(
"Response" = c(log(MMR.NMS.Data$MMR)), "Fit" = c(predict(MMR.m3)),
"Mass" = c(log(MMR.NMS.Data$Mass))
)
}
MMR.NMS.Fit.m1 <- MMR.NMS.Fit.m1 %>%
mutate(model = "Model 1")
MMR.NMS.Fit.m2 <- MMR.NMS.Fit.m2 %>%
mutate(model = "Model 2")
MMR.NMS.Fit.m3 <- MMR.NMS.Fit.m3 %>%
mutate(model = "Model 3")
MMR.NMS.Fit.All <- rbind(MMR.NMS.Fit.m1, MMR.NMS.Fit.m2, MMR.NMS.Fit.m3)
MMR.NMS.Fits <- ggplot(data = MMR.NMS.Fit.All, aes(
x = Fit, y = Response, fill = model, colour = model,
linetype = model
)) +
theme_bw() +
my.theme +
geom_point(pch = 21, size = 3, alpha = 0.5) +
geom_smooth(method = "lm") +
xlab("Predicted log Resting Metabolic Rate \n (log mg O2/h)") +
ylab("Measured log Resting \n Metabolic Rate \n (log mg O2/h)") +
scale_fill_manual(values = c(
"mediumseagreen",
"mediumorchid", "cornflowerblue"
)) +
scale_colour_manual(values = c(
"black",
"black", "black"
)) +
theme(legend.title = element_blank())
print(MMR.NMS.Fits)
ggsave("log Peak Metabolic Rate Model Fit Comparisons.jpeg", dpi = 1000, limitsize = TRUE)
# Very subtle differences between models.Amalgamating mass-independant results.
require("plyr")
R2.Dataframe <- ldply(lapply(list(CTM.model.m1, RMR.m1, RMR.m2, RMR.m3, MMR.m1, MMR.m2, MMR.m3), r.squaredGLMM), data.frame)
{
CTM.NMS.DF1 <- subset(CTM.Models, model == "Model 1") %>%
mutate("df" = 8) %>%
mutate("Model Number" = "Model 1") %>%
mutate("Marginal R2" = R2.Dataframe[1, 1]) %>%
mutate("Conditional R2" = R2.Dataframe[1, 2]) %>%
mutate("AIC" = summary(CTM.model.m1)$AIC)
CTM.NMS.DF1$Response <- "Critical Thermal Maximum"
RMR.NMS.DF1 <- subset(RMR.NMS.Models, model == "Model 1") %>%
mutate("df" = 8) %>%
mutate("Model Number" = "Model 1") %>%
mutate("Marginal R2" = R2.Dataframe[2, 1]) %>%
mutate("Conditional R2" = R2.Dataframe[2, 2]) %>%
mutate("AIC" = summary(RMR.m1)$AICtab[1])
RMR.NMS.DF2 <- subset(RMR.NMS.Models, model == "Model 2") %>%
mutate("df" = 8) %>%
mutate("Model Number" = "Model 2") %>%
mutate("Marginal R2" = R2.Dataframe[3, 1]) %>%
mutate("Conditional R2" = R2.Dataframe[3, 2]) %>%
mutate("AIC" = summary(RMR.m2)$AICtab[1])
RMR.NMS.DF3 <- subset(RMR.NMS.Models, model == "Model 3") %>%
mutate("df" = 8) %>%
mutate("Model Number" = "Model 3") %>%
mutate("Marginal R2" = R2.Dataframe[4, 1]) %>%
mutate("Conditional R2" = R2.Dataframe[4, 2]) %>%
mutate("AIC" = summary(RMR.m3)$AICtab[1])
RMR.NMS.DF <- rbind(RMR.NMS.DF1, RMR.NMS.DF2, RMR.NMS.DF3)
RMR.NMS.DF$Response <- "Resting Metabolic Rate"
RMR.NMS.DF$Response <- "Resting Metabolic Rate"
RMR.NMS.DF$Response <- "Resting Metabolic Rate"
MMR.NMS.DF1 <- subset(MMR.NMS.Models, model == "Model 1") %>%
mutate("df" = 8) %>%
mutate("Model Number" = "Model 1") %>%
mutate("Marginal R2" = R2.Dataframe[5, 1]) %>%
mutate("Conditional R2" = R2.Dataframe[5, 2]) %>%
mutate("AIC" = summary(MMR.m1)$AICtab[1])
MMR.NMS.DF2 <- subset(MMR.NMS.Models, model == "Model 2") %>%
mutate("df" = 8) %>%
mutate("Model Number" = "Model 2") %>%
mutate("Marginal R2" = R2.Dataframe[6, 1]) %>%
mutate("Conditional R2" = R2.Dataframe[6, 2]) %>%
mutate("AIC" = summary(MMR.m2)$AICtab[1])
MMR.NMS.DF3 <- subset(MMR.NMS.Models, model == "Model 3") %>%
mutate("df" = 8) %>%
mutate("Model Number" = "Model 3") %>%
mutate("Marginal R2" = R2.Dataframe[7, 1]) %>%
mutate("Conditional R2" = R2.Dataframe[7, 2]) %>%
mutate("AIC" = summary(MMR.m3)$AICtab[1])
MMR.NMS.DF <- rbind(MMR.NMS.DF1, MMR.NMS.DF2, MMR.NMS.DF3)
MMR.NMS.DF$Response <- "Peak Metabolic Rate"
MMR.NMS.DF$Response <- "Peak Metabolic Rate"
MMR.NMS.DF$Response <- "Peak Metabolic Rate"
}
AllData.NMS <- rbind(CTM.NMS.DF1, RMR.NMS.DF, MMR.NMS.DF)
AllData.NMS <- AllData.NMS %>%
select("Response", "Model Number", "term", "estimate", "std.error", "statistic", "df", "p.values", "AIC", "Marginal R2", "Conditional R2")
colnames(AllData.NMS) <- c("Response Variable", "Model Number", "Model Parameter", "Estimate (beta)", "Std.Error", "T Statistic", "DF", "p", "AIC", "Marginal R2", "Conditional R2")
New.Parameters <- c()
for (i in 1:nrow(AllData.NMS)) {
New.Parameters[i] <- gsub("\n", "", AllData.NMS$"Model Parameter"[i])
}
AllData.NMS$"Model Parameter" <- New.Parameters
rownames(AllData.NMS) <- c(1:nrow(AllData.NMS))
View(AllData.NMS)
write.csv(AllData.NMS, "Mass Independant Model Output Results - Final.csv")Producing final plots
## Critical thermal maximum
# Top model
CTM.model.m1 <- lme(CTM ~ Offspring.Acc, random = list(~ 1 | Mother.ID, ~ 1 | Father.ID), weights = varIdent(~ Mass | Father.Acc * Offspring.Acc), data = CTM.NMS.Data)
# Means and plots
require("emmeans")
Mat_mmeans <- as.data.frame(emmeans(CTM.model.m1, ~Offspring.Acc, data = CTM.NMS.Data, nesting = NULL, nesting.order = FALSE, CIs = TRUE, plot = FALSE))
Mat_expand <- expand.grid("Mother.Acc" = c(1, 2), "Offspring.Acc" = c(1, 2))
Mat_mmeans <- Mat_expand %>%
mutate(
"yvar" = ifelse(Offspring.Acc == 1, Mat_mmeans$emmean[1], Mat_mmeans$emmean[2]),
"SE" = ifelse(Offspring.Acc == 1, Mat_mmeans$SE[1], Mat_mmeans$SE[2])
)
Pat_mmeans <- as.data.frame(emmeans(CTM.model.m1, ~Offspring.Acc, data = CTM.NMS.Data, nesting = NULL, nesting.order = FALSE, CIs = TRUE, plot = FALSE))
Pat_expand <- expand.grid("Father.Acc" = c(1, 2), "Offspring.Acc" = c(1, 2))
Pat_mmeans <- Pat_expand %>%
mutate(
"yvar" = ifelse(Offspring.Acc == 1, Pat_mmeans$emmean[1], Pat_mmeans$emmean[2]),
"SE" = ifelse(Offspring.Acc == 1, Pat_mmeans$SE[1], Pat_mmeans$SE[2])
)
Mat_dw <- ggplot(Mat_mmeans, aes(x = Mother.Acc, y = yvar, fill = factor(Offspring.Acc))) +
geom_jitter(
geom = "point", data = CTM.NMS.Data, mapping = aes(x = Mother.Acc, y = CTM), size = 1,
pch = 21, alpha = 0.3, position = position_dodge(width = 0.5)
) +
geom_errorbar(
mapping = aes(x = Mother.Acc, ymin = yvar - SE * 1.96, ymax = yvar + SE * 1.96),
size = 1, width = 0.2, colour = "black", position = position_dodge(width = 0.5)
) +
geom_point(size = 4, colour = "black", pch = 21, position = position_dodge(width = 0.5)) +
scale_x_discrete(labels = c("Cold Mothers", "Warm Mothers")) +
ylab("Critical Thermal Maximum (°C)") +
scale_fill_manual(
values = c("dodgerblue3", "magenta"),
name = "",
labels = c(" Cold Offspring", "Warm Offspring")
) +
my.theme +
theme_bw() +
theme(axis.title.x = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/CTM_mmeans_mothers.jpg", Mat_dw,
width = 6, height = 6, dpi = 800
)
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/CTM_mmeans_mothers.pdf", Mat_dw,
width = 6, height = 6, dpi = 800
)
Pat_dw <- ggplot(Pat_mmeans, aes(x = Father.Acc, y = yvar, fill = factor(Offspring.Acc))) +
geom_jitter(
geom = "point", data = CTM.NMS.Data, mapping = aes(x = Father.Acc, y = CTM), size = 1,
pch = 21, alpha = 0.3, position = position_dodge(width = 0.5)
) +
geom_errorbar(
mapping = aes(x = Father.Acc, ymin = yvar - SE * 1.96, ymax = yvar + SE * 1.96),
size = 1, width = 0.2, colour = "black", position = position_dodge(width = 0.5)
) +
geom_point(size = 4, colour = "black", pch = 21, position = position_dodge(width = 0.5)) +
scale_x_discrete(labels = c("Cold Fathers", "Warm Fathers")) +
ylab("Critical Thermal Maximum (°C)") +
scale_fill_manual(
values = c("dodgerblue3", "magenta"),
name = "",
labels = c("Cold Offspring", "Warm Offspring")
) +
my.theme +
theme_bw() +
theme(axis.title.x = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/CTM_mmeans_fathers.jpg", Pat_dw,
width = 6, height = 6, dpi = 800
)
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/CTM_mmeans_fathers.pdf", Pat_dw,
width = 6, height = 6, dpi = 800
)
# Residual plots
CTM_Resid <- lme(CTM ~ Mass, random = list(~ 1 | Mother.ID, ~ 1 | Father.ID), weights = varIdent(~ Mass | Father.Acc * Offspring.Acc), data = CTM.NMS.Data)
CTM.NMS.Data$Resid_CTM <- residuals(CTM_Resid, type = "response")
Mat_resid_CTM <- ggplot(CTM.NMS.Data, aes(x = Mother.Acc, y = Resid_CTM, fill = Offspring.Acc)) +
geom_point(size = 1, pch = 21, alpha = 0.3, position = position_dodge(width = 0.5)) +
stat_summary(
geom = "errorbar", fun.data = "mean_cl_boot", size = 1,
width = 0.2, colour = "black", position = position_dodge(width = 0.5)
) +
stat_summary(
geom = "point", fun = "mean", size = 4, colour = "black",
pch = 21, position = position_dodge(width = 0.5)
) +
scale_x_discrete(labels = c("Cold Mothers", "Warm Mothers")) +
ylab("Critical Thermal Maximum (Residual °C)") +
scale_fill_manual(
values = c("dodgerblue3", "magenta"),
name = "",
labels = c("Cold Offspring", "Warm Offspring")
) +
my.theme +
theme_bw() +
theme(axis.title.x = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/CTM_residuals_mothers.jpg", Mat_resid_CTM,
width = 6, height = 6, dpi = 800
)
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/CTM_residuals_mothers.pdf", Mat_resid_CTM,
width = 6, height = 6, dpi = 800
)
Pat_resid_CTM <- ggplot(CTM.NMS.Data, aes(x = Father.Acc, y = Resid_CTM, fill = Offspring.Acc)) +
geom_point(size = 1, pch = 21, alpha = 0.3, position = position_dodge(width = 0.5)) +
stat_summary(
geom = "errorbar", fun.data = "mean_cl_boot", size = 1,
width = 0.2, colour = "black", position = position_dodge(width = 0.5)
) +
stat_summary(
geom = "point", fun = "mean", size = 4, colour = "black",
pch = 21, position = position_dodge(width = 0.5)
) +
scale_x_discrete(labels = c("Cold Fathers", "Warm Fathers")) +
ylab("Critical Thermal Maximum (Residual °C)") +
scale_fill_manual(
values = c("dodgerblue3", "magenta"),
name = "",
labels = c("Cold Offspring", "Warm Offspring")
) +
my.theme +
theme_bw() +
theme(axis.title.x = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/CTM_residuals_fathers.jpg", Pat_resid_CTM,
width = 6, height = 6, dpi = 800
)
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/CTM_residuals_fathers.pdf", Pat_resid_CTM,
width = 6, height = 6, dpi = 800
)
## Resting MO2
# Top model
RMR.m1 <- lmer(log(RMR) ~ log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID), REML = FALSE, data = RMR.NMS.Data)
Mat_mmeans <- as.data.frame(emmeans(RMR.m1, ~Mother.Acc, data = RMR.NMS.Data, nesting = NULL, nesting.order = FALSE, CIs = TRUE, plot = FALSE))
Mat_expand <- expand.grid("Mother.Acc" = c(1, 2), "Offspring.Acc" = c(1, 2))
Mat_mmeans <- Mat_expand %>%
mutate(
"yvar" = ifelse(Mother.Acc == 1, exp(Mat_mmeans$emmean[1]),
exp(Mat_mmeans$emmean[2])
),
"UCL" = ifelse(Mother.Acc == 1, exp(Mat_mmeans$emmean[1] + 1.96 * Mat_mmeans$SE[1]),
exp(Mat_mmeans$emmean[2] + 1.96 * Mat_mmeans$SE[2])
),
"LCL" = ifelse(Mother.Acc == 1, exp(Mat_mmeans$emmean[1] - 1.96 * Mat_mmeans$SE[1]),
exp(Mat_mmeans$emmean[2] - 1.96 * Mat_mmeans$SE[2])
),
)
RMR_mmeans1 <- as.data.frame(emmeans(RMR.m1, ~Mass, cov.reduce = FALSE, data = RMR.NMS.Data, nesting = NULL, nesting.order = FALSE, CIs = TRUE, plot = FALSE))
Pat_expand <- expand.grid("Father.Acc" = c(1, 2), "Offspring.Acc" = c(1, 2))
Pat_mmean <- Pat_expand %>%
mutate(
"yvar" = exp(RMR_mmeans1$emmean),
"UCL" = exp(RMR_mmeans1$emmean + 1.96 * RMR_mmeans1$SE),
"LCL" = exp(RMR_mmeans1$emmean - 1.96 * RMR_mmeans1$SE)
)
RMR_father_plot <- ggplot(Pat_mmean, aes(x = Father.Acc, y = yvar, fill = factor(Offspring.Acc))) +
geom_jitter(
geom = "point", data = RMR.NMS.Data, mapping = aes(x = Father.Acc, y = RMR),
size = 1, pch = 21, alpha = 0.3, position = position_dodge(width = 0.5)
) +
geom_errorbar(
mapping = aes(x = Father.Acc, ymin = LCL, ymax = UCL),
size = 1, width = 0.2, colour = "black", position = position_dodge(width = 0.5)
) +
geom_point(size = 4, colour = "black", pch = 21, position = position_dodge(width = 0.5)) +
scale_x_discrete(labels = c("Cold Fathers", "Warm Fathers")) +
ylab(expression(Routine ~ MO[2] ~ (mg ~ O[2] ~ h^{
"-1"
}))) +
scale_fill_manual(
values = c("dodgerblue3", "magenta"),
name = "",
labels = c("Cold Offspring", "Warm Offspring")
) +
my.theme +
theme_bw() +
theme(axis.title.x = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/RMR1_mmeans_by_fathers.jpg", RMR_father_plot, width = 6, height = 6, dpi = 800)
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/RMR1_mmeans_by_fathers.pdf", RMR_father_plot, width = 6, height = 6, dpi = 800)
RMR_mother_plot <- ggplot(Mat_mmeans, aes(x = Mother.Acc, y = yvar, fill = factor(Offspring.Acc))) +
geom_jitter(
geom = "point", data = RMR.NMS.Data, mapping = aes(x = Mother.Acc, y = RMR),
size = 1, pch = 21, alpha = 0.3, position = position_dodge(width = 0.5)
) +
geom_errorbar(
mapping = aes(x = Mother.Acc, ymin = LCL, ymax = UCL),
size = 1, width = 0.2, colour = "black", position = position_dodge(width = 0.5)
) +
geom_point(size = 4, colour = "black", pch = 21, position = position_dodge(width = 0.5)) +
scale_x_discrete(labels = c("Cold Mothers", "Warm Mothers")) +
ylab(expression(Routine ~ MO[2] ~ (mg ~ O[2] ~ h^{
"-1"
}))) +
scale_fill_manual(
values = c("dodgerblue3", "magenta"),
name = "",
labels = c("Cold Offspring", "Warm Offspring")
) +
my.theme +
theme_bw() +
theme(axis.title.x = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/RMR1_mmeans_by_mothers.jpg", RMR_mother_plot, width = 6, height = 6, dpi = 800)
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/RMR1_mmeans_by_mothers.pdf", RMR_mother_plot, width = 6, height = 6, dpi = 800)
# Residual plots
RMR_Resid <- lmer(log(RMR) ~ log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID), REML = FALSE, data = RMR.NMS.Data)
RMR.NMS.Data$Resid_RMR <- residuals(RMR_Resid, type = "response")
Mat_resid_RMR <- ggplot(RMR.NMS.Data, aes(x = Mother.Acc, y = Resid_RMR, fill = Offspring.Acc)) +
geom_point(size = 1, pch = 21, alpha = 0.3, position = position_dodge(width = 0.5)) +
stat_summary(
geom = "errorbar", fun.data = "mean_cl_boot", size = 1,
width = 0.2, colour = "black", position = position_dodge(width = 0.5)
) +
stat_summary(
geom = "point", fun = "mean", size = 4, colour = "black",
pch = 21, position = position_dodge(width = 0.5)
) +
scale_x_discrete(labels = c("Cold Mothers", "Warm Mothers")) +
ylab(expression(~Residual ~ Routine ~ MO[2] ~ (Log ~ mg ~ O[2] ~ h^{
"-1"
}))) +
scale_fill_manual(
values = c("dodgerblue3", "magenta"),
name = "",
labels = c("Cold Offspring", "Warm Offspring")
) +
geom_hline(yintercept = 0, colour = "black", size = 0.5) +
my.theme +
theme_bw() +
theme(axis.title.x = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/RMR_residuals_mothers_log.jpg", Mat_resid_RMR, width = 6, height = 6, dpi = 800)
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/RMR_residuals_mothers_log.pdf", Mat_resid_RMR, width = 6, height = 6, dpi = 800)
Pat_resid_RMR <- ggplot(RMR.NMS.Data, aes(x = Father.Acc, y = Resid_RMR, fill = Offspring.Acc)) +
geom_point(size = 1, pch = 21, alpha = 0.3, position = position_dodge(width = 0.5)) +
stat_summary(
geom = "errorbar", fun.data = "mean_cl_boot", size = 1,
width = 0.2, colour = "black", position = position_dodge(width = 0.5)
) +
stat_summary(
geom = "point", fun = "mean", size = 4, colour = "black",
pch = 21, position = position_dodge(width = 0.5)
) +
scale_x_discrete(labels = c("Cold Fathers", "Warm Fathers")) +
ylab(expression(~Residual ~ Routine ~ MO[2] ~ (Log ~ mg ~ O[2] ~ h^{
"-1"
}))) +
scale_fill_manual(
values = c("dodgerblue3", "magenta"),
name = "",
labels = c("Cold Offspring", "Warm Offspring")
) +
geom_hline(yintercept = 0, colour = "black", size = 0.5) +
my.theme +
theme_bw() +
theme(axis.title.x = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/RMR_residuals_fathers_log.jpg", Pat_resid_RMR, width = 6, height = 6, dpi = 800)
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/RMR_residuals_fathers_log.pdf", Pat_resid_RMR, width = 6, height = 6, dpi = 800)
# Peak MO2
MMR.m1 <- lmer(log(MMR) ~ log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID), REML = FALSE, data = MMR.NMS.Data)
# Means and plots
# Plotting marginal means
Mat_mmeans <- as.data.frame(emmeans(MMR.m1, ~Offspring.Acc, data = MMR.NMS.Data, nesting = NULL, nesting.order = FALSE, CIs = TRUE, plot = FALSE))
Mat_expand <- expand.grid("Mother.Acc" = c(1, 2), "Offspring.Acc" = c(1, 2))
Mat_mmeans <- Mat_expand %>%
mutate(
"yvar" = ifelse(Offspring.Acc == 1, exp(Mat_mmeans$emmean[1]),
exp(Mat_mmeans$emmean[2])
),
"UCL" = ifelse(Offspring.Acc == 1, exp(Mat_mmeans$emmean[1] + 1.96 * Mat_mmeans$SE[1]),
exp(Mat_mmeans$emmean[2] + 1.96 * Mat_mmeans$SE[2])
),
"LCL" = ifelse(Offspring.Acc == 1, exp(Mat_mmeans$emmean[1] - 1.96 * Mat_mmeans$SE[1]),
exp(Mat_mmeans$emmean[2] - 1.96 * Mat_mmeans$SE[2])
)
)
Pat_mmeans <- as.data.frame(emmeans(MMR.m1, ~Offspring.Acc, data = MMR.NMS.Data, nesting = NULL, nesting.order = FALSE, CIs = TRUE, plot = FALSE))
Pat_expand <- expand.grid("Father.Acc" = c(1, 2), "Offspring.Acc" = c(1, 2))
Pat_mmeans <- Pat_expand %>%
mutate(
"yvar" = ifelse(Offspring.Acc == 1, exp(Pat_mmeans$emmean[1]),
exp(Pat_mmeans$emmean[2])
),
"UCL" = ifelse(Offspring.Acc == 1, exp(Pat_mmeans$emmean[1] + 1.96 * Pat_mmeans$SE[1]),
exp(Pat_mmeans$emmean[2] + 1.96 * Pat_mmeans$SE[2])
),
"LCL" = ifelse(Offspring.Acc == 1, exp(Pat_mmeans$emmean[1] - 1.96 * Pat_mmeans$SE[1]),
exp(Pat_mmeans$emmean[2] - 1.96 * Pat_mmeans$SE[2])
)
)
Mat_MMR <- ggplot(Mat_mmeans, aes(x = Mother.Acc, y = yvar, fill = factor(Offspring.Acc))) +
geom_jitter(
geom = "point", data = MMR.NMS.Data, mapping = aes(x = Mother.Acc, y = MMR), size = 1,
pch = 21, alpha = 0.3, position = position_dodge(width = 0.5)
) +
geom_errorbar(
mapping = aes(x = Mother.Acc, ymin = LCL, ymax = UCL),
size = 1, width = 0.2, colour = "black", position = position_dodge(width = 0.5)
) +
geom_point(size = 4, colour = "black", pch = 21, position = position_dodge(width = 0.5)) +
scale_x_discrete(labels = c("Cold Mothers", "Warm Mothers")) +
ylab(expression(Peak ~ MO[2] ~ (mg ~ O[2] ~ h^{
"-1"
}))) +
scale_fill_manual(
values = c("dodgerblue3", "magenta"),
name = "",
labels = c(" Cold Offspring", "Warm Offspring")
) +
my.theme +
theme_bw() +
theme(
axis.title.x = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/MMR_mmeans_mothers.jpg", Mat_MMR,
width = 6, height = 6, dpi = 800
)
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/MMR_mmeans_mothers.pdf", Mat_MMR,
width = 6, height = 6, dpi = 800
)
Pat_MMR <- ggplot(Pat_mmeans, aes(x = Father.Acc, y = yvar, fill = factor(Offspring.Acc))) +
geom_jitter(
geom = "point", data = MMR.NMS.Data, mapping = aes(x = Father.Acc, y = MMR), size = 1,
pch = 21, alpha = 0.3, position = position_dodge(width = 0.5)
) +
geom_errorbar(
mapping = aes(x = Father.Acc, ymin = LCL, ymax = UCL),
size = 1, width = 0.2, colour = "black", position = position_dodge(width = 0.5)
) +
geom_point(size = 4, colour = "black", pch = 21, position = position_dodge(width = 0.5)) +
scale_x_discrete(labels = c("Cold Fathers", "Warm Fathers")) +
ylab(expression(Peak ~ MO[2] ~ (mg ~ O[2] ~ h^{
"-1"
}))) +
scale_fill_manual(
values = c("dodgerblue3", "magenta"),
name = "",
labels = c("Cold Offspring", "Warm Offspring")
) +
my.theme +
theme_bw() +
theme(axis.title.x = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/MMR_mmeans_fathers.jpg", Pat_MMR,
width = 6, height = 6, dpi = 800
)
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/MMR_mmeans_fathers.pdf", Pat_MMR,
width = 6, height = 6, dpi = 800
)
# Residual plots
MMR_Resid <- lmer(log(MMR) ~ log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID), REML = FALSE, data = MMR.NMS.Data)
MMR.NMS.Data$Resid_MMR <- residuals(MMR_Resid, type = "response")
Mat_resid_MMR <- ggplot(MMR.NMS.Data, aes(x = Mother.Acc, y = Resid_MMR, fill = Offspring.Acc)) +
geom_point(size = 1, pch = 21, alpha = 0.3, position = position_dodge(width = 0.5)) +
stat_summary(
geom = "errorbar", fun.data = "mean_cl_boot", size = 1,
width = 0.2, colour = "black", position = position_dodge(width = 0.5)
) +
stat_summary(
geom = "point", fun = "mean", size = 4, colour = "black",
pch = 21, position = position_dodge(width = 0.5)
) +
scale_x_discrete(labels = c("Cold Mothers", "Warm Mothers")) +
ylab(expression(Peak ~ MO[2] ~ (Log ~ Residual ~ mg ~ O[2] ~ h^{
"-1"
}))) +
scale_fill_manual(
values = c("dodgerblue3", "magenta"),
name = "",
labels = c("Cold Offspring", "Warm Offspring")
) +
geom_hline(yintercept = 0, colour = "black", size = 0.5) +
my.theme +
theme_bw() +
theme(axis.title.x = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/MMR_residuals_mothers_Log.jpg", Mat_resid_MMR, width = 6, height = 6, dpi = 800)
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/MMR_residuals_mothers_Log.pdf", Mat_resid_MMR,
width = 6, height = 6, dpi = 800
)
Pat_resid_MMR <- ggplot(MMR.NMS.Data, aes(x = Father.Acc, y = Resid_MMR, fill = Offspring.Acc)) +
geom_point(size = 1, pch = 21, alpha = 0.3, position = position_dodge(width = 0.5)) +
stat_summary(
geom = "errorbar", fun.data = "mean_cl_boot", size = 1,
width = 0.2, colour = "black", position = position_dodge(width = 0.5)
) +
stat_summary(
geom = "point", fun = "mean", size = 4, colour = "black",
pch = 21, position = position_dodge(width = 0.5)
) +
scale_x_discrete(labels = c("Cold Fathers", "Warm Fathers")) +
ylab(expression(Peak ~ MO[2] ~ (Log ~ Residual ~ mg ~ O[2] ~ h^{
"-1"
}))) +
scale_fill_manual(
values = c("dodgerblue3", "magenta"),
name = "",
labels = c("Cold Offspring", "Warm Offspring")
) +
geom_hline(yintercept = 0, colour = "black", size = 0.5) +
my.theme +
theme_bw() +
theme(axis.title.x = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/MMR_residuals_fathers.jpg", Pat_resid_MMR, width = 6, height = 6, dpi = 800)
ggsave("/home/joshk/Documents/Trout/Brook_Trout_Revised_Figures/MMR_residuals_fathers.pdf", Pat_resid_MMR, width = 6, height = 6, dpi = 800)